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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

The  dynamic  response  of  a  submersible  vehicle  operating  at  the  extremes 
of  its  operational  envelope  is  becoming  increasingly  important  in  order  to  en- 
hance vehicle  operations.  Traditionally,  dynamic  stability  of  motion  is  stud- 
ied using  eigenvalue  analysis  where  the  equations  of  motion  are  linearized 
around  nominal  straight  line  level  flight  paths  (Arentzen  &:  Mandel,  1960), 
(Clayton  &;  Bishop,  1982),  (Feldman,  1987).  Under  certain  simplified  as- 
sumptions, a  simple  criterion  Gv  >  0  can  be  obtained  where  the  stability 
index  Gv  is  function  of  the  hydrodynamic  coefficients  in  heave  and  pitch. 
Values  for  the  stability  index  can  be  computed  by, 

Mw{Zq+m) 
Gv  =  1~        ZwMq        •  (1) 

This  index  is  analogous  to  the  familiar  stability  coefficient  for  horizontal 
plane  maneuvering  (Lewis,  1989)  and  can  be  thought  of  as  a  high  speed  ap- 
proximation where  the  effect  of  the  metacentric  restoring  moment  is  minimal. 
If  the  value  of  Gv  is  greater  than  zero,  the  vehicle  is  dynamically  stable.  As 
we  point  out  in  the  next  chapter  though,  this  is  only  a  sufficient,  and  rather 
conservative  condition  for  stability.  It  is  not  a  necessary  condition  in  the 
sense  that  Gv  <  0  indicates  instability  at  infinite  forward  speed.  It  is  quite 
possible  that  at  normal  operating  speeds  the  vehicle  might  be  directionally 
stable.   Furthermore,  Gv  <  0  indicates  a  divergent  loss  of  stability  which  is 


quite  uncommon  in  the  vertical  plane.  Most  modern  submarines  exhibit  a 
flutter-like  instability  at  high  speed,  which  can  not  be  analyzed  using  the 
above  simplified  index.  Divergent  motions  may  develop  in  combined  six 
degrees  of  freedom  (Papoulias  et  al  1993)  and  their  occurrence  can  not  be 
analyzed  by  a  single  stability  index. 


B.  OBJECTIVES  AND  OUTLINE 

In  this  work  we  examine  the  problem  of  stability  of  motion  with  controls 
fixed  in  the  vertical  plane,  with  particular  emphasis  on  the  mechanism  of 
loss  of  stability  of  straight  line  motion.  The  surge  equation  is  decoupled 
from  heave/pitch  through  a  perturbation  series  approach  (Bender  Sz  Orszag, 
1978).  It  is  shown  that  loss  of  stability  occurs  in  the  form  of  generic  bifurca- 
tions to  periodic  solutions  (Guckenheimer  &  Holmes,  1983).  Taylor  expan- 
sions and  center  manifold  approximations  are  employed  in  order  to  isolate 
the  main  nonlinear  terms  that  influence  system  response  after  the  initial  loss 
of  stability  (Hassard  &:  Wan,  1978).  Integral  averaging  is  performed  in  order 
to  combine  the  nonlinear  terms  into  a  design  stability  coefficient  (Chow  & 
Mallet-Paret,  1977).  Special  attention  is  paid  to  the  study  of  the  quadratic 
drag  terms  as  they  constitute  some  of  the  main  nonlinear  terms  of  the  equa- 
tions of  motion.  This  is  a  unique  feature  of  the  open  loop  problem.  In  similar 
studies  of  loss  of  stability  under  closed  loop  depth  control  (Bateman,  1993) 
it  was  found  that  the  main  damping  mechanism  is  provided  through  the  use 


of  control  susrafces.  The  difficulty  associated  with  the  nonsmoothness  of  the 
absolute  value  nonlinearities  of  the  quadratic  drag  forces  is  dealt  with  by 
employing  the  concept  of  generalized  gradient  (Clarke,  1983).  This  has  the 
advantage  of  keeping  the  linear  terms  constant,  unlike  the  linear/cubic  ap- 
proximation typically  used  in  ship  roll  motion  studies  (Dalzell,  1978),  where 
the  linear  damping  coefficient  is  a  function  of  the  assumed  amplitude  of  mo- 
tion. 

Vehicle  modeling  in  this  work  follows  standard  notation  (Gertler  &:  Hagen, 
1976),  (Smith  et  al  1978),  and  numerical  results  are  presented  for  the  DARPA 
SUBOFF  model  (Roddy,  1990)  for  which  a  set  of  hydrodynamic  coefficients 
and  geometric  properties  is  available.  Furthermore,  the  baseline  vehicle  is 
marginally  stable  with  controls  fixed  under  normal  operating  conditions  and 
can  serve  as  a  prime  example  for  the  techniques  described  in  this  work.  The 
model  has  been  experimentally  validated  for  angles  of  attack  on  the  hull 
between  ±15  deg.,  while  the  constant  coefficient  approximation  introduces 
very  little  error  in  time  domain  simulations  (Tinker,  1978).  Unless  otherwise 
mentioned,  all  results  in  this  work  are  presented  in  standard  dimensionless 
form  with  respect  to  the  vehicle  length  I  =  4.26  m,  and  nominal  forward 
speed  U  =  2.44  m/sec.  All  angular  deflections  are  shown  in  degrees. 


n.  PROBLEM  FORMULATION 


A.  EQUATIONS  OF  MOTION 

Assuming  that  vehicle  motion  is  restricted  in  the  vertical  plane,  the  math- 
ematical model  consists  of  the  coupled  nonlinear  heave  and  pitch  equations 
of  motion.  In  a  moving  coordinate  frame  fixed  at  the  vehicle's  geometrical 
center,  Newton's  equations  of  motion  for  a  port/starboard  symmetric  and 
neutrally  buoyant  vehicle  are  expressed  in  dimensionless  form  as  follows, 

m(w  -uq  -  zGq7  -  xGq)  =  Zqq  +  Z^w  +  Zqq  +  Zww 

/•nose 

-CD  \        b(x){w  -  xq)\w  -  xq\  dx  +  Zs8  ,  (2) 

./tail 

Iyq  +  mzG(u  +  wq)  —  mxG(w  —  uq)  =  Mqq  +  M^w  +  Mqq  +  Mww 

/•nose 

+Cd  I        b(x)(w  —  xq)\w  —  xq\x  dx 

./tail 

-xGBW cos  0  -  zGBW sin  d  +  Ms6  ,  (3) 

where  xGb  =  xg  —  xb,  zGb  =  zG  —  zb  ,  and  the  rest  of  the  symbols  are  based 
on  standard  notation  and  are  explained  in  Table  1.  Without  loss  of  generality 
we  can  assume  that  zb  =  xb  =  0,  so  that  xGb  =  xG  and  zGb  =  zG.  The  cross 
flow  integral  terms  in  these  equations  become  very  important  for  high  angles 
of  attack  maneuvering,  where  they  provide  the  primary  motion  damping. 
The  drag  coefficient,  Cd,  is  assumed  to  be  constant  throughout  the  vehicle 
length  for  simplicity.    This  does  not  affect  the  qualitative  properties  of  the 


results  that  follow.  The  vehicle  pitch  rate  is, 

9  =  q.  (4) 

Dynamic  coupling  between  surge  and  heave/pitch  is  present  due  to  coordi- 
nate coupling  as  a  result  of  the  nonzero  metacentric  height.  Therefore,  pitch 
and  heave  motions  must  be  studied  together  with  surge, 

mu  +  mwq  -  mxGq2  +  rnzGq  =  Xqqq2  +  X^u  -f  Xwqwq  +  Xwww2 

+Xuuu2  +  Xnnn2  +  X6S62  ,  (5) 

where  we  assume  that  both  resistance  and  propulsive  forces  are  proportional 
to  the  square  of  the  speed  or  the  propeller  revolutions,  respectively. 

In  analyzing  controls  fixed  stability  of  motion,  the  case  8  =  0  is  examined 
first.  The  steady  state  solutions  of  the  equations  are  determined  by  w  =  q  = 
u  =  $  =  q0  =  0,  where  subscript  0  indicates  variable  value  at  steady  state. 
Substituting  these  conditions  in  (2)  we  get, 

Zww0  -  CoAvWolwol  =  0  ,  (6) 

where, 

/•nose 

K  =  /        b(x)  dx  ,  (7) 

./tail 

is  the  "waterplane"  area.  Since  Zw  <  0,  equation  (7)  admits  only  one  solu- 
tion, namely  w0  =  0.  Equation  (3)  then  yields, 

tan0o  =  -~,  (8) 

zgb 

while  (5)  is  used  to  determine  the  nominal  forward  speed,  Uq. 


TABLE  1:  NOMENCLATURE 


a 

dummy  independent  variable 

oo 

steady  state  value  of  a 

<*ij 

expansion  coefficients  of  23  in  terms  of  Z\ ,  22 

b(x) 

local  beam  of  the  hull 

cD 

quadratic  drag  coefficient 

7 

regularization  parameter 

S 

stern  plane  deflection 

£ 

perturbation  parameter,  e  =  (m4) 

€ 

criticality  difference,  e  =  u  —  uc 

4 

vehicle  mass  moment  of  inertia 

K 

nonlinear  stability  coefficient 

A 

system  eigenvalue 

m 

vehicle  mass 

M 

pitch  moment 

M0 

derivative  of  M  with  respect  to  a 

9 

pitch  rate 

(*,*) 

polar  coordinates  of  Z\ ,  z2 

T 

transformation  matrix  of  x  to  z 

e 

pitch  angle 

u 

forward  speed 

Ue 

critical  value  of  u 

w 

heave  velocity 

X 

state  variables  vector,  x  =  [9,w,q] 

X 

surge  force 

xa 

derivative  of  X  with  respect  to  a 

(*b,zb) 

body  fixed  coordinates  of  vehicle  center  of  buoyancy 

{*g,zg) 

body  fixed  coordinates  of  vehicle  center  of  gravity 

xgb 

center  of  gravity /center  of  buoyancy  separation,  xq  —  %b 

ZGB 

vehicle  metacentric  height,  Zq  —  z& 

z 

state  variables  vector  in  its  normal  form 

«1,  22 

critical  coordinates  of  z 

^3 

stable  coordinate  of  z 

z 

heave  force 

za 

derivative  of  Z  with  respect  to  a 

Wo 

imaginary  part  of  critical  pair  of  eigenvalues 

B.  REDUCTION  OF  ORDER 

The  linearized  surge,  heave,  and  pitch  equations  of  motion  in  the  vicinity 
of  the  nominal  point  are, 

(m  -  X„)u  +  mzGq  =  2Xuuu  ,  (9) 

(m  -  Zv)w  -  (mxG  +  Zq)  =  (Zq  +  m)q  +  Zww  ,  (10) 

(Iy  -  Mq)q  +  mzGu  -  {mxG  +  M^w  =  Mww  +  (Mq  -  mxG)q 

+M9e,  (ii) 

where, 

Mb  =  xGB  W  sin  0O  -  zgb  W  cos  60  ,  (12) 

is  the  hydrostatic  restoring  moment  coefficient.  The  characteristic  equation 
of  (4),  (9),  (10),  and  (11)  is  obtained  as, 

(_AXA  +  BX)(A2\Z  +  B2X2  +  C2A  +  D3)  +  A3X4  +  B3\3  =  0  ,  (13) 

where, 

Ax  =  m  —  X„  , 

B\  =  2XUU  , 

A2  =  (m  -  Zv)(Iv  -  Mq)  -  (Zq-  +  ma^XM,*  +  miG)  , 

52  =  -Zw(/y  -  Mq)  -  (m  -  Z^)(Mq  -  mxG)  -  (Zq  +  m){M^  +  ma;G) 

-M^Z,,  +mxG)  , 

C2  =  -M«(m  -  Zw)  +  ZW(M,  -  mxG)  -  MW(Z,  +  m)  , 

D7  =  MeZw  , 


A3     =     {mzG)2{m  -  Zy,)  , 
B3     =    —(raze)  Zw  . 

It  can  be  seen  that  the  parameter  (tuzq)  is  responsible  for  surge  and  heave, 

pitch  coupling.  For  zq  =  0,  equation  (13)  decouples  into  the  surge  eigenvalue 

A  =  B\/A\  and  the  classical  cubic  characteristic  equation  for  the  vertical 

plane.     It  should  be  mentioned  that  the  effect  of  the  forward  speed  u  is 

embedded  into  the  definition  for  the  dimensionless  vehicle  weight  W  through, 

W 
^pu2L2 

If  we  introduce  a  smallness  parameter, 

e  =  K)2  ,  (15) 

we  can  rewrite  (13)  in  the  form, 

(A  +  ea)\4  +  (B  +  e/?)A3  +  CX2  +  DX  +  E  =  0  ,  (16) 

where  in  terms  of  previously  defined  coefficients, 

A  =  -AxOt.  > 

B  =  -A1B7+B1A2  , 

C  =  -AXC7+BXB2  , 

D  =  -AxD2  +  BlC2  , 

E  =  BXD2  , 

a  =  m  —  Z„  , 

0  =  -Zw  . 
8 


Following  (Bender  Sz  Orszag,  1978)  we  expand  the  solutions  of  (16)  in  a 
regular  perturbation  series, 

\  =  \0  +  \le  +  O{e2)  ,  (17) 

where  A0  is  an  eigenvalue  for  e  =  0  (uncoupled  surge  or  heave /pitch),  Xx  is 
the  first  order  correction  due  to  dynamic  coupling,  and  0(e2)  contains  second 
and  higher  order  terms  in  e.  Substituting  (17)  into  (16)  we  can  get, 


Ag(aAo+/3) 
AAX30  +W\20  +  2CA0  +  D 


a = a0 -  AAX3^::r:^  ,  ^ + o^) .       m 


It  can  be  seen  that  the  correction  term  is  very  small  compared  to  the  un- 
coupled root  as  evidenced  by  the  Ajj  term.  This  is  particularly  true  when 
A0  nears  zero;  i.e.,  close  to  a  bifurcation  point.  Therefore,  loss  of  stability 
can  be  studied  by  analyzing  the  heave/pitch  equations  decoupled  from  surge. 
The  characteristic  equation  then  becomes, 


A2X3  +  B7\2  +  C2A  +  D2  =  0  .  (19) 

Plots  of  the  system  eigenvalues  at  nominal  speed  versus  zq  are  shown  in 
Figures  1  through  4.  The  surge  eigenvalues  is  real  and  negative  throughout 
the  range  of  zq,  while  the  heave/pitch  eigenvalues  are  real  for  small  values 
of  zq.  The  two  larger  real  heave/pitch  eigenvalues  coalesce  into  a  complex 
conjugate  pair  whose  real  part  crosses  zero  for  a  certain  value  of  zq.  Within 
the  accuracy  of  Figures  1  through  3,  the  eigenvalues  A  are  identical  to  those 
computed  by  either  the  coupled  or  the  uncoupled  system,  or  the  perturbation 
equations  (18).    There  is  a  very  small  difference  for  the  surge  eigenvalue  as 


shown  in  Figure  4,  but  again  the  agreement  between  coupled  and  uncoupled 
systems  is  very  good. 


C.  DEGREE  OF  STABILITY 

The  eigenvalues  of  the  reduced  order  system  (18)  designate  stability  or 
instability  of  motion.  A  single  measure  of  stability,  the  "degree  of  stability", 
e„,  can  be  defined  as  the  maximum  real  part  of  all  three  eigenvalues  of  (18). 
This  measures  the  slowest  exponential  convergence  to  the  equilibrium  when 
negative  and  the  fastest  exponential  divergence  from  the  equilibrium  when 
positive.  For  all  numerical  calculations  in  this  work,  the  degree  of  stability 
corresponds  to  the  real  part  of  a  complex  conjugate  pair  of  eigenvalues.  Typ- 
ical results  are  presented  in  Figures  5  through  7.  Figure  5  shows  the  degree 
of  stability  versus  xqb  constant  dimensionless  speed  u0  =  0.5  and  different 
values  of  zqq.  Similar  results  are  shown  in  Figure  6  for  constant  zqb  =  0.015 
and  different  values  of  uq.  Finally,  a  three  dimensional  representation  is 
depicted  in  Figure  7. 


D.  CRITICAL  SPEED 

The  parameter  value  where  the  real  part  of  the  complex  conjugate  pair  of 
eigenvalues  shown  in  the  previous  figures  crosses  zero  defines  the  point  where 
linear  stability  is  lost.    This  critical  point  can  be  computed  by  considering 
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Figure  1:  System  eigenvalues  versus  zq 
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Figure  2:  Complex  conjugate  heave/pitch  eigenvalues  versus  zq 
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Figure  3:  Real  heave/pitch  eigenvalues  versus  zq 
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Figure  4:  Surge  eigenvalue  versus  zq 
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Figure  5:  Degree  of  stability  versus  xqb  for  uq  =  0.5  and  different  values  of  zqb 
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Figure  6:  Degree  of  stability  versus  xGB  for  zqB  -  0.015  and  different  values  of  uq 
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Figure  7:  Degree  of  stability  versus  xqb  and  uo  for  zqb  =  0.015 
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Figure  8:  Critical  speed  u,.  versus  xq  for  different  values  of  zqb 
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Figure  9:  Critical  speed  Ue  versus  xq  and  zq 
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equation  (19).  Routh's  criterion  applied  to  this  cubic  yields  A2D2  =  BiC2 
which  can  be  solved  for  the  dimensionless  weight, 

■A-2  -i/2,1   —  -"2^2,1 

where, 

C7fi     =     Zw(Mq  -  mxG)  -  Mw(Zq  +  m)  , 
C2il     =     (m  -  Z„){zGB  cos  0O  -  %gb  sin  60)  > 

■^2,1      =      ZW(XGB  Sin  ^0  —  2TGB  cos  #o)  • 

The  value  of  the  critical  speed  ue  can  then  be  evaluated  from  (20)  and  (14). 
Typical  results  are  presented  in  Figures  8  and  9,  nondimensionalized  with 
respect  to  nominal  vehicle  speed  2.44  m/sec  and  length  4.26  m.  Vertical  plane 
motions  are  stable  for  forward  speeds  less  than  the  critical  speed.  It  can  be 
seen  that  stability  is  increasing  with  increasing  zq  while  xq  =  0  is  the  most 
conservative  condition  for  stability.  Therefore,  a  vehicle  which  is  stable  when 
properly  trimmed  will  remain  stable  for  off-trim  conditions.  For  comparison, 
we  note  that  the  simple  stability  coefficient  Gv,  defined  in  equation  (1),  is 
monotonically  decreasing  and  becomes  more  negative  for  decreasing  xG,  as 
shown  in  Figure  10.  Thus  it  would  have  predicted  unstable  motions  for  the 
entire  range  of  parameters  shown  in  Figures  8  and  9.  For  completeness,  the 
value  of  the  steady  state  pitch  angle,  0o  (in  degrees),  is  shown  in  Figure  11 
versus  xGb  using  zGb  as  the  parameter.  This  pitch  angle  is  computed  from 
equation  (8). 
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Figure  10:  Stability  coefficient  Gv  versus  xgb 
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III.  BIFURCATION  ANALYSIS 


A.  INTRODUCTION 

In  all  cases  of  stability  loss  of  the  previous  chapter,  one  pair  of  com- 
plex conjugate  eigenvalues  of  the  corresponding  eigenvalue  problem  crosses 
transversally  the  imaginary  axis.  A  situation  like  this  in  which  a  certain 
parameter  is  varied  such  that  the  real  part  of  one  pair  of  complex  conjugate 
eigenvalues  of  the  linearized  system  matrix  crosses  zero,  results  in  the  sys- 
tem leaving  its  steady  state  in  an  oscillatory  manner.  This  loss  of  stability 
is  called  Hopf  bifurcation  and  generically  occurs  in  either  supercritical  or 
subcritical  form.  In  the  supercritical  case,  stable  limit  cycles  are  generated 
after  the  nominal  straight  line  motion  loses  its  stability.  The  amplitudes  of 
these  limit  cycles  are  continuously  increasing  as  the  parameter  distance  from 
its  critical  value  is  increased.  For  small  values  of  this  criticality  distance  the 
resulting  limit  cycle  is  of  small  amplitude  and  differs  little  from  the  initial 
nominal  state.  In  the  subcritical  case,  however,  stable  limit  cycles  are  gener- 
ated before  the  nominal  state  loses  its  stability.  Therefore,  depending  on  the 
initial  conditions  it  is  possible  to  diverge  away  from  the  nominal  straight  line 
path  and  converge  towards  a  limit  cycle  even  before  the  nominal  motion  loses 
its  stability.  This  means  that  in  the  subcritical  Hopf  bifurcation  case  the  do- 
main of  attraction  of  the  nominal  state  is  decreasing  and  in  fact  it  shrinks 
to  zero  as  the  critical  point  is  approached.  Random  external  disturbances  of 
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sufficient  magnitude  can  throw  the  vehicle  off  to  an  oscillatory  steady  state 
even  though  the  nominal  state  may  still  remain  stable.  After  the  nominal 
state  becomes  unstable,  a  discontinuous  increase  in  the  magnitude  of  motions 
is  observed  as  there  exist  no  simple  stable  nearby  attractors  for  the  vehicle 
to  converge  to.  Distinction  between  these  two  qualitatively  different  types  of 
bifurcation  is,  therefore,  essential  in  design.  The  computational  procedure 
requires  higher  order  approximations  in  the  equations  of  motion  and  is  the 
subject  of  this  chapter. 


B.  THIRD  ORDER  EXPANSIONS 

The  nonlinear  heave/pitch  equations  of  motion  (2),  (3),  and  (4)  are  writ- 
ten in  the  form, 

8    =    <?,  (21) 

w     =    anw  +  ax2q  +  a,i3(xGB  cos  0  +  zGB  sin  0)  +  (^(w^q) 

+**(»,«),  (22) 

q    =    a21w  +  022q  +  a-x&GB  cos  9  +  zGB  sin  0)  +  dq(w,  q) 

+c2(w,q)  ,  (23) 


where 


Dv     =    (m  -  Z*)(Iy  -  Mq)  -  [mxG  +  Zq)(mxG  +  M«)  , 
auDv     =    (Iy  -  Mq)Zw  +  (mxG  +  Zq)Mw  , 
al2Dv     =    (Iy  -  Mq)(m  +  Zq)  +  (mxG  +  Zq)(Mq  -  mxG)  , 
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a13Dv     =     -{mxG  +  Zj)W  , 
d2\Dv     =    (m  -  Z^)MW  +  (mxG  +  M^)ZW  , 
«22^v     =    (wi  -  Z^M,  -  mxG)  +  (mxc  +  M^)(m  +  Z,)  , 
<l23-Ov     =     —  (wi  -  Z^)^  , 
c^(tx;,g)Dv     =    (7y  -  M^  +  (mxG  +  Z^)/,  , 
dq(w,q)Dv     =    (m  -  Z^)/,  +  (miG  +  M„)IW  , 
Oi(i0,g).Dv     =     (Jv  -  M^)mzGq7  -  (mxG  +  Zq)mzGwq  , 
c2(w,q)Dv     =     — (m  -  Z^)mzGwq  +  (mxG  +  M^mzcq7  , 
and  /„,,  i,  are  the  cross  flow  integrals 

/«,     =    Cd  /        6(a?)(tw  -  xg)|tu  -  xg|<Zs  ,  (24) 

•/tail 

/•nose 

■^     =    C*£>  /        b(x)(w  —  xq)\w  —  xglxdsc  .  (25) 

•/tail 

The  system  of  equations  (21)  through  (23)  is  written  in  the  compact  form 

x  =  Ax+g(x),  (26) 

where 

x  =  [0,«;ig],  (27) 

is  the  three  state  variables  vector,  and  A  is  the  linearized  sytem  matrix  eval- 
uated at  the  nominal  point  Xq.  The  term  g(x)  contains  all  nonlinear  terms 
of  the  equations.  Hopf  bifurcation  analysis  can  be  performed  by  isolating 
the  primary  nonlinear  terms  in  g(x).  Keeping  terms  up  to  third  order,  we 
can  write 

g(x)  =  g(2)(x)  +  g(3)(x).  (28) 
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Using  equations  (21)  through  (25),  the  various  terms  in  (28)  can  be  written 


as, 


9?     =    0, 


9?     =    {ly  ~  Mi)mzGq7  -  (mxG  +  Z4)mzGwq  +  ^{w^q)  ,  (29) 

§3       =     ~(m  -  Zu,)inzGwq  +  (mxG  +  M„}mzGq*  +  ^{vj,q)  , 


and 


9?     =    0, 

02°     =    <&3)K  9)  +  iai3(*GB  sin  90  -  zGB  cos  eo)03  ,  (30) 

03       =    tf\w,  ^)  +  \<h*{?GB  sin  #o  -  ^GB  cos  00)82  • 

Expansion  in  Taylor  series  of  <£„,,  <fq  requires  expansion  of  the  cross  flow 
integrals  Iw,  Iq,  which  require  the  Taylor  series  of 

/K)  =  Ml  •  (3i) 

This  expression  can  be  converted  into  an  analytic  function  using  Dalzell's 
approximation  (Dalzell,  1978), 

which  is  derived  by  a  least  squares  fit  of  an  odd  series  over  some  assumed 
range  of  £,  namely  —  £c  <  £  <  £..  This  approximation,  which  is  shown  in 
Figure  12,  has  been  extensively  used  in  ship  roll  motion  studies  and  is  very 
useful  for  its  intended  purpose.  However,  in  the  present  problem  it  suffers 
from  several  major  drawbacks: 
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•  It  introduces  a  linear  term  which  depends  on  the  assumed  range  of 
motion,  and  it  renders  the  critical  speed  function  of  the  vehicle  motions. 

•  The  cubic  term,  which  is  ultimately  responsible  for  the  Hopf  bifurcation 
analysis,  is  a  function  of  the  assumed  range  of  vehicle  motions  which 
can  not  be  known  in  advance. 

•  As  Figure  12  demonstrates,  the  slope  of  the  actual  curve  at  the  origin 
is  significanly  different  than  the  approximation,  which  would  make  the 
bifurcation  results  unreliable. 

Instead  of  Dalzell's  approximation,  we  employ  the  concept  of  generalized 
gradient  (Clarke,  1983),  which  is  used  in  the  study  of  control  systems  in- 
volving discontinuous  or  non-smooth  functions.  In  this  way  we  approximate 
the  gradient  of  a  non-smooth  function  at  a  discontinuity  by  a  map  equal  to 
the  convex  closure  of  the  limiting  gradients  near  the  discontinuity.  In  our 
problem  we  write, 

/«)  =  &I&I  +  2|6|«  -  fc)  +  sign(£o)(£  -  £o)2  +  /(3)(0  ,  (33) 

as  the  Taylor  series  epansion  of  /(£)  near  £0-  The  sign  function  in  (33)  can 
be  approximated  by, 

sign(£0)  =  limtanh  (  — J   .  (34) 

A  graphical  representation  of  the  approximation  (34)  is  shown  in  Figure 
13.    The  quantity  7  is  a  small  regularization  parameter  and  is  used  in  the 
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Figure  12:   Graphical  representation  of  DalzelTs  approximation  of  £|£|  versus  £/£. 
Solid  curve  is  the  exact  expression  and  dotted  curve  is  the  approximation  (32) 
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next  section  for  the  proper  normalization  of  the  results.   Using  (34),  we  can 
approximate  /(£)  in  the  vicinity  of  £0  =  0  by, 

Ml  *  ^3  •  (35) 

Since 

£  (—►  w  —  xq  ,  (36) 

we  can  express  the  non-smooth  cross  flow  integral  terms  by, 

Iw     =     ~-{E0w3  -  3E,w7q  +  3E7wq7  -  E*q3)  ,  (37) 

67 

/,     =     ^(Etw3  -3E7w7q  +  3E3wq7  -  EAq3)  ,  (38) 

67 

where 

/>no«e 

Ei  =  /        xlb{x)dx  ,  (39) 

./tail 

are  the  moments  of  the  vehicle  "waterplane"  area. 


C.  COORDINATE  TRANSFORMATIONS 

Using  the  previous  second  and  third  order  Taylor  series  expansions,  equa- 
tion (26)  is  written  in  the  form, 

x  =  Ax  +  g(2>(x)+g(3)(x)  .  (40) 

If  T  is  the  matrix  of  eigenvectors  of  A  evaluated  at  the  critical  point  u  =  uc, 
the  linear  change  of  coordinates, 

x  =  Tz  ,     z  =  T_1x  ,  (41) 
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Figure  13:  Graphical  representation  of  the  sign  function  and  its  hyperbolic  tangent 
approximation  (34).  Solid  curve  corresponds  to  7  =  0.1  and  dotted  curve  corre- 
sponds to  7  =  0.01 
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transforms  system  (40)  into  its  normal  coordinate  form, 

z  =  T_1  ATz  +  T_1g(2)(Tz)  +  T_1g(3)(Tz)  .  (42) 

At  the  Hopf  bifurcation  point,  matrix  T_1AT  takes  the  form, 

T'AT  = 

where  a>0  is  the  imaginary  part  of  the  critical  pair  of  eigenvalues,  and  the 
remaining  eigenvalue  p  is  negative.  For  values  of  u  close  to  the  bifurcation 
poit  ttc,  matrix  T-IAT  becomes, 


0 

-u0 

0  " 

u>0 

0 

0 

0 

0 

P  . 

T_1AT  = 


a'e  -(wo  +  w'e)         0 

(u;0  +  w'e)  a'e  0 

0  0  p  +  j/c 


where  e  denotes  the  criticality  difference 


e  =  u  —  uf 


(43) 


and 


a     = 


U)       = 


V     = 


derivative  of  the  real  part  of  the  critical  eigenvalue 

with  respect  to  e  , 

derivative  of  the  imaginary  part  of  the  critical  eigenvalue 

with  respect  to  e  , 

derivative  of  p  with  respect  to  e  . 
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D.  CENTER  MANIFOLD  EXPANSIONS 

Due  to  continuity,  the  eigevalue  p  +  p'e  remains  negative  for  small  nonzero 
values  of  e.  Therefore,  the  coordinate  z3  corresponds  to  a  negative  eigenvalue 
and  is  asymptotically  stable.  Center  manifold  theory  predicts  that  the  rela- 
tionship between  the  critical  coordinates  z1?  z2  and  the  stable  coordinate  z3 
is  at  least  of  quadratic  order.  We  can  then  write  z3  as, 


Z3   =  aU2j    -|-  G-\lZ\Zl  +  OL77Z7    , 


(44) 


where  the  coefficients,  a^-,  in  the  quadratic  center  manifold  expansion  (44) 
need  to  be  determined.  By  differentiating  equation  (44)  we  obtain, 


i3  =  2a11z1z1  +  0:12(21 22  +  ziz2)  +  2a22z2z2  • 


(45) 


We  substitute  i\  =  —u)0z2  and  z2  —  wo^i  from  equation  (42)  into  (45),  and 
we  obtain 


i3  =  al2uj0z\  +  2(a22  -  "liVo^i^  —  clX2WqZ17 


(46) 


The  third  equation  of  (42)  is  written  as, 


i3=^  +  [T-'g«(Tz)]w), 


(47) 


where  terms  up  to  second  order  have  been  kept.   If  we  denote  the  elements 
of  T  and  T_1  by, 

T  =  [my]  ,     T"1  =  [ny]  ,  (48) 


then 


T-ig(2)(Tz)  = 


d2 
d$ 
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where 

dx     =    nl7{£nzl  +  i76ziz7  +  £77z27)  +  nl3(t3Szl  +  l^zxz7  +  £37z77)  ,    (49) 

d2      =      n77(£7bzl  +  ll*ZXZ7  +  ^27^2)  +  n23  (4s  z\  +  t^ZxZ-i.  +  ^37^2)  >     (50) 
<4       =      n32  (4b  2?   +  4s  ZX  Z2  +  ^7  222  )  +  "33  (4»5  z\   +  lM  Zj  Z2   +  ^7  z\  )   .      (51 ) 

Expressions  for  the  coefficients  4j  axe  given  in  the  Fortran  programs  in  the 
appendix. 

Equation  (47)  then  becomes 

z3  =  pz3  -f  d$  ,  (52) 

and  substituting  (44)  and  (51)  into  (52)  we  get, 

i3     =    (pan  +  n32^25  +  "33^35)^1  +  {v0L\7  +  "32^26  +  ntt£zG)z\Zi 

+(P«22  +  7132^27  +  n3347)-Z2    '  (53) 

Comparing  coefficients  of  (46)  and  (53)  we  get, 

-p<XU  +  <*>0<*12      =     "32^5  +  Tl3345  ,  (54) 

— 2u>0aii  —  pa\7  +  2u>oQ!22     =    ^32^26  +  ™33^36  ,  (55) 

— U;0«i2  —  PCC22       =      "32^27  +  ™33^37   •  (56) 

Solution  of  the  system  of  linear  equations  (54)  through  (56)  yields  the  coef- 
ficients in  the  center  manifold  expansion  (44). 
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E.  AVERAGING 

Using  the  previous  Taylor  expansions  and  center  manifold  approxima- 
tions, we  can  write  the  reduced  two-dimensional  system  that  describes  the 
center  manifold  flow  of  (42)  in  the  form, 

zx     =    a'ezl  -  (u>0  +  u/e)z2  +  ^1(^1,22)  1  (57) 

z2     =    (uq  +  oj'e)zl  +  a'ez7  +  F2(zx,z2)  ,  (58) 


where 


^i(*i»«a)     =    rii^i  +r122122  +r13ziz7  +  ruz7 

+Pn*i  +  P12Z1Z2  +  Pi34  ,  (59) 

F2{zuz2)     =    r^+mzlzt+r^zl+rnz} 

+P7i*i  +P21Z1Z7  +P23Z7  ?  (60) 


and 


r, 


ij  =  ni7£2j  +  ru3£3j  ,     i  =  1, 2  ,     j  =  1, . . . ,  4  ,  (61) 

A;  =  ni24fc  +  ",-34*  ,     t  =  l,2,     ;  =  1,2,3  ,      fc  =  j  +  4.     (62) 

If  we  introduce  polar  coordinates  in  the  form, 

zx  =  R  cos  <f>  ,      z2  =  i2  sin  <£  ,  (63) 

we  can  use  (57)  and  (58)  to  produce  an  equation  describing  the  rate  of  change 
of  the  radial  coordinate  R, 

R  =  a'eR  +  P(<f>)R3  +  Q{<f>)R2  .  (64) 
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This  equation  contains  one  variable,  R,  which  is  slowly  varying  in  time,  and 
another  variable,  (j>,  which  is  a  fast  variable.  Therefore,  equation  (64)  can  be 
averaged  over  one  complete  cycle  in  (f)  to  produce  an  equation  with  constant 
coefficients  and  similar  stability  properties, 

R  =  aeR  +  KR3  +  LR2  ,  (65) 

where 

K    =     J-  [^  P(<f>)d<f>=l(3ru+rl3+r77+3r74),  (66) 

Z7T  JO 

L    =    ^Jo*rQ(<f>)d<j>  =  0.  (67) 

Therefore,  the  averaged  equation  (65)  becomes 

R  =  a'eR  +  KR3  .  (68) 


F.  LIMIT  CYCLE  ANALYSIS 

Equation  (68)  admits  two  steady  state  solutions,  one  at  R  =  0  which 
corresponds  to  the  trivial  equilibrium  solution  at  zero,  and  one  at 


R«  =  \J-K€'  (69) 

This  equilibrium  solution  corresponds  to  a  periodic  solution  or  limit  cycle  in 
the  cartesian  coordinates  Zi,  z2.  For  this  limit  cycle  to  exist,  the  quantity  Ro 
must  be  a  real  number.  In  our  case  a'  is  always  positive,  since  the  system 
loses  its  stability;  i.e.,  the  real  part  of  the  critical  pair  of  eigenvalues  changes 
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from  negative  to  positive,  for  increasing  u.     Therefore,  existence  of  these 
periodic  solutions  depends  on  the  value  of  K.  Specifically, 

•  if  K  <  0,  periodic  solutions  exist  for  e  >  0  or  u  >  tic,  and 

•  if  K  >  0,  periodic  solutions  exist  for  e  <  0  or  u  <  uc. 

The  characteristic  root  of  (68)  in  the  vicinity  of  (69)  is 

/3  =  -2a'e  ,  (70) 

and  we  can  see  that 

•  if  periodic  solutions  exist  for  u  >  uc  they  are  stable,  and 

•  if  periodic  solutions  exist  for  u  <  uc  they  are  unstable. 

The  period  of  these  periodic  solutions  can  be  estimated  as  follows.  Equa- 
tions (57),  (58),  and  (63)  produce  an  equation  in  <f>  similar  to  (64), 

j>  =  u;0  +  u/e  +  F{<f>)R2  +  G{4>)R  .  (71) 

The  averaged  form  of  (71)  is 

<£  =  o;0  +  Je  +  MB?  ,  (72) 

where 

M=±r  F(4>)  d<f>  =  \(Zrn  +  r73  -  r12  -  3r14)  .  (73) 

*/0 


The  limit  cycle  period  can  be  computed  by  substituting  (69)  into  (72), 


2tt  2tt  /        JK  -  a'M 


u;o  +  u/e  +  MKq       ojq  \  ujqK 

36 


e)  +  O  (C2)   .         (74) 


G.  RESULTS  AND  DISCUSSION 

Typical  results  of  the  nonlinear  stability  coefficient  K  are  shown  in  Fig- 
ures 14  and  15.  Figure  14  presents  a  plot  of  K  •  7  versus  xq  for  zq  =  0.015 
and  for  different  values  of  Cd-  It  should  be  emphasized  that  the  use  of  K  -7 
is  more  meaningful  than  the  use  of  K,  since  it  properly  accounts  for  the 
use  of  the  regularization  parameter  7  as  seen  from  equations  (35)  and  (68). 
Numerical  evidence  demonstrates  that  all  curves  K  •  7  versus  xq  converge 
for  7  — ►  0.  For  practical  purposes,  values  of  7  smaller  than  0.001  produce 
identical  results.  The  results  of  Figure  14  demonstrate  the  profound  effect 
that  the  quadratic  drag  coefficient  Cd  has  on  stability  of  limit  cycles.  All 
Hopf  bifurcations  are  supercritical  (K  <  0),  and  they  become  stronger  su- 
percritical as  Cd  is  increased.  It  is  worth  noting  that  results  for  Cd  =  0 
produce  subcritical  behavior,  K  >  0,  which  is  clearly  incorrect.  Thus,  ne- 
glecting the  effects  of  Cd  would  have  produce  entirely  wrong  results  in  the 
present  problem.  Figure  15  shows  a  plot  of  K  •  7  versus  xq  for  Cd  =  0.5 
and  different  values  of  the  metacentric  height  zq.  It  can  be  seen  that,  the 
bifurcations  become  stronger  supercritical  as  initial  stability  zq  is  increased. 

The  bifurcation  analysis  results  are  verified  by  direct  numerical  simula- 
tions shown  in  Figures  16  through  18.  Figure  16  shows  the  results  of  two 
numerical  simulations  for  two  values  of  nominal  speed  tto  in  terms  of  the 
vehicle  pitch  angle  9  (in  degrees)  versus  time  (in  seconds).  The  critical  value 
of  speed,  uc,  is  about  0.495  as  can  be  seen  from  Figure  8,  while  the  other 
parameters  in  the  simulations  were  Cd  =  0.5,  zq  =  0.015,  and  xq  =  0.    It 
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Figure  14:  K  -j  versus  xG  for  zq  =  0.015  and  different  values  of  the  drag  coefficient 
CD 
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Figure  15:    K  •  7  versus  xq  for  Co  =  0.5  and  different  values  of  the  metacentric 
height  zq 
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Figure  16:  Time  histories  (9,t)  for  CD  =  0.5,  zq  =  0.015,  xG  =  0,  and  two  different 
values  of  nominal  speed 
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Figure  17:  Time  histories  (0,  t)  for  Co  —  0.5,  zq  —  0.015,  xq  —  0,  and  a  range  of 
nominal  forward  speeds 
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Figure  18:  Limit  cycle  amplitudes  from  the  results  of  Figure  17 
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can  be  seen  that  convergence  to  zero  is  ensured  for  Uq  <  uc  and  convergence 
to  a  limit  cycle  occurs  for  uq  >  uc.  This  indicates  supercritical  behavior  as 
shown  before.  A  selection  of  time  histories  is  shown  in  Figure  17  for  a  range 
of  forward  speeds  and  the  same  parameters  as  in  Figure  16.  The  same  initial 
disturbance,  0  =  5  degrees,  was  introduced  at  t  =  0  for  all  simulations.  It 
can  be  seen  that  the  amplitude  of  limit  cycles  increases  as  the  distance  of  u 
from  uc  is  increasing.  The  rate  of  convergence  of  solutions  to  their  limit  cy- 
cles is  also  increasing,  while  their  period  remains  essentially  constant.  These 
results  are  summarized  in  Figure  18,  where  the  amplitudes  of  the  numeri- 
cally computed  limit  cycles  are  plotted  versus  u$.  The  behavior  is  clearly 
supercritical,  which  agrees  with  our  findings  of  the  bifurcation  analysis. 
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IV.  BIAS  EFFECTS 


A.  LOSS  OF  STABILITY 

Stability  analysis  of  motions  at  a  nonzero  angle  of  attack  can  be  per- 
formed by  first  introducing  some  bias  into  the  steady  state  solution  and  its 
perturbations.  This  can  be  achieved  by  maintaining  a  nonzero  dive  plane 
angle  at  nominal.  In  this  case  the  steady  state  solutions  are  <fo  =  0,  and  itfo, 
0O  are  computed  from, 

ZwwQ  -  CD E0w0 \w0\  +  Zs6  =  0  ,  (75) 

MwwQ  +  CdExwq\w0\  -  W(xGB  cos  0O  +  zgb  sin  d0)  +  M68  =  0  .  (76) 

The  coefficients  E0,  E\  are  computed  from  (39).    In  order  to  solve  (75)  we 
observe  that  when  8  >  0,  then  w0  <  0, 


-Zw  -  JZI-4CDE0ZS6 
Wo  =  ?r   F '  (77) 


and  when  8  <  0,  then  w0  >  0, 


-Zw  -  JZ*+4CDEoZs8 
™o  =  r    = •  (78) 

The  angle  90  can  then  be  computed  from  equation  (76). 

Linearization  of  the  equations  of  motion  in  the  neighborhood  of  the  above 
equilibrium  point  produces  the  linear  system, 

(m  -  Z^)w  -  (mxG  +  Zj)q  —  (Zw  -  2CDE0\w0\)w 
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+{Zq+m  +  2CDEl\w0\)qi  (79) 

-(M*  +  mxG)w  +  (/„  -  M,j)9  =  {Mw  +  2CD£!  |w0|)u; 

+(M,  -  mxG  -  mzGw0  -  2CDE7\w0\)q 

+W(xGB  sin  O0  -  zgb  cos  6o)0  ,  (80) 

$  =  q,  (81) 

where  the  variables  w,  6  are  understood  as  small  deviations  from  their  equi- 
librium values.  Numerical  solution  of  the  generalized  eigenvalue  problem  (79) 
through  (81)  yields  the  critical  speed  values  where  the  nominal  equilibrium 
solution  becomes  unstable. 


B.  ANALYSIS  OF  HOPF  BIFURCATIONS 

It  can  be  numerically  verified  that  the  above  calculations  for  the  new 
critical  speed  result  in  a  loss  of  stability  in  the  form  of  Hopf  bifurcations,  as 
for  the  6  =  0  case.  These  Hopf  bifurcations  can  be  analyzed  using  the  same 
general  methodology  that  was  developed  in  the  previous  chapter.  The  non- 
linear expansions  are  like  equations  (21)  to  (23)  with  the  following  changes: 
The  substitutions 

Zw     — ►     Zw  —  2CdEqWq  , 
Mw    -*    Mw  +  2CDElw0  , 
Zq  +m    -*     Zq  +m  +  2CDEiW0  , 
Mq  —  mxG     — ►     Mq  —  mxG  —  mzGw0  —  2CdE2W0  , 
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axe  assumed  in  the  definition  of  coefficients  a^.  Furthermore,  equation  (33) 
prduces  2nd  order  contributions  due  to  w0  ^  0,  as  well  as  3rd  order.  Using 
(33)  we  can  compute  the  2nd  order  expansions  of  d$  and  S2^  in  (29)  using, 

4,2)     =    CD(E0w2  -2Elwq  +  E7q2),  (82) 

/J2)     =    CD{Exw2  -  2E7wq  +  E3q2)  .  (83) 

These  equations  are  valid  for  wq  >  0  or  8  <  0.  For  8  >  0  the  signs  of  E{ 
must  be  switched.  The  third  order  expansions  i^3)  and  J^3^  are  the  same 
as  in  equations  (37)  and  (38).  Using  these  additional  terms,  the  nonlinear 
stability  coefficient  K  can  be  computed  in  the  same  way  as  in  the  previous 
chapter. 


C.  RESULTS  AND  DISCUSSION 

Typical  results  for  8  ^  0  are  shown  in  Figures  19  through  24.  Figure 
19  shows  the  equilibrium  pitch  angle  60  (in  degrees)  versus  xq  for  different 
values  of  8  from  —5  to  5  degrees  with  increments  of  1  degree.  Solid  curves 
correspond  to  positive  8  and  dashed  curves  to  negative.  Figure  20  shows  the 
degree  of  stability  for  the  equilibrium  points  of  Figure  19,  while  Figure  21 
presents  the  degree  of  stability  in  a  three  dimensional  view.  It  can  be  seen 
that  positive  and  negative  values  of  8  have  almost  identical  stability  charac- 
teristics. Furthermore,  the  degree  of  stability  becomes  more  negative  as  the 
absolute  value  of  8  is  increased,  which  means  that  we  expect  a  wider  domain 
of  stability  in  this  case.  This  is  verified  by  the  critical  speed  plots  shown  in 
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Figure  19:  Pitch  angle  9Q  versus  xG  for  zq  =  0.015,  u<,  =  0.5,  and  different  values 
of* 
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Figure  20:    Degree  of  stability  versus  xq  for  zq  —  0.015,  uo  =  0.5,  and  different 
values  of  6 
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Figure  21:  Degree  of  stability  versus  xq  and  uo  for  zq  —  0.015  and  two  values  of  8 
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Figure  22:  Critical  speed  Uc  versus  xq  for  zq  =0.015  and  different  values  of  8 
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Figure  23:  Critical  speed  u«.  versus  xq  and  6  for  zq  —  0.015 
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Figure  24:   Nonlinear  stability  coefficient  K  versus  xq  for  Co  =  0.5,  zq  =  0.015, 
and  different  values  of  6 
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Figures  22  and  23.  The  critical  speed  is  minimum  for  8  =  0  and  it  increases 
monotonically  with  increasing  absolute  value  of  8.  This  stabilizing  effect  of 
asymmetry  (bias)  remains  approximately  true  for  the  nonlinear  analysis,  as 
demonstrated  by  the  results  of  Figure  24.  It  can  be  seen  that  the  nonlinear 
stability  coefficient  K  becomes  more  negative  as  8  is  decreasing  from  zero. 
For  increasing  8,  K  becomes  less  negative  but  the  difference  from  the  6  =  0 
calculations  appears  to  be  very  small.  Therefore,  limit  cycle  stability  is  not 
significantly  affected  by  the  bias  effects  that  are  induced  by  small  nonzero 
dive  plane  angles. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  presented  a  comprehensive  nonlinear  study  of  straight  line  stabil- 
ity of  motion  of  submersibles  in  the  dive  plane  under  open  loop  conditions. 
A  systematic  perturbation  analysis  demonstrated  that  the  effects  of  surge  in 
heave/pitch  are  small  and  can  be  neglected.  Primary  loss  of  stability  was 
shown  to  occur  in  the  form  of  Hopf  bifurcations  to  periodic  solutions.  The 
critical  speed  were  instability  occurs  was  computed  in  terms  of  metacentric 
height,  longitudinal  separation  of  the  centers  of  buoyancy  and  gravity,  and 
the  dive  plane  angle.  Analysis  of  the  periodic  solutions  that  resulted  from  the 
Hopf  bifurcations  was  accomplished  through  Taylor  expansions,  up  to  third 
order,  of  the  equations  of  motion.  A  consistent  approximation,  utilizing  the 
generalized  gradient,  was  used  to  study  the  non-analytic  quadratic  cross  flow 
integral  drag  terms.  The  results  indicated  that  loss  of  stability  occurs  always 
in  the  form  of  supercritical  Hopf  bifurcations  with  stable  limit  cycles.  It  was 
shown  that  this  is  mainly  due  to  the  stabilizing  effect  of  the  drag  forces  at 
high  angles  of  attack.  As  a  recommendation  for  further  research,  we  suggest 
a  nonlinear  analysis  of  coupled  motions  in  six  degrees  of  freedom. 
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APPENDIX 


The  following  is  a  list  and  description  of  the  computer  programs  used  in  this 
thesis.  The  programs  are  written  in  FORTRAN  or  MATLAB.  Complete 
printouts  of  the  programs  follow  after  the  list. 

•  PERTURB.M 

MATLAB  program  for  performing  surge/heave/pitch  perturbation  anal- 
ysis. 

•  CRITO.M 

MATLAB  program  for  calculating  the  critical  speed  for  6  =  0. 

•  CRITJDELTA.M 

MATLAB  program  for  calculating  the  critical  speed  for  6^0. 

•  HOPF_0.FOR 

FORTRAN  program  for  calculating  the  nonlinear  stability  coefficient 
K  for  6  =  0. 

•  HOPF.DELTA.FOR 

FORTRAN  program  for  calculating  the  nonlinear  stability  coefficient 
K  for  8  ^  0. 

•  SIM.FOR 

FORTRAN  program  for  simulating  the  equations  of  motion. 
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7,  Program  perturb. m 

'/,  This  file  applies  concepts  of  perturbation  analysis  theory, 
'/,  in  order  to  prove  that  the  3  by  3  system  describing  motion 
'/,  in  the  vertical  plane,  decouples  from  the  4  by  4  one. 


rho  = 

1. 

94; 

8   = 

32 

•  2; 

L 

13 

.9792; 

ndi  = 

0. 

5*rho*L"2 

nd2  = 

0. 

5*rho*L~3 

nd3  = 

0. 

5*rho*L~4 

nd4  = 

0. 

5*rho*L~5 

iy  = 

561.32/nd4; 

zqdot 

= 

-6.33e-4; 

zvdot 

= 

-1.4529e-2; 

zq 

= 

7.545e-3; 

zv 

= 

-1.391e-2; 

mqdot 

= 

-8.8e-4; 

mvdot 

= 

-5.61e-4; 

mq 

= 

-3.702e-3; 

mw 

= 

1.0324e-2; 

cd 

= 

0.015; 

zb 

= 

C/L; 

m 

= 

1556. 2363/ (g*nd2); 

xudot 

= 

-0.05*m; 

xb 

= 

0/L; 

uo 

= 

4; 

w 

= 

1556 . 2363 . / (ndl . *uo . ~2) ; 

b 

= 

w; 

xg    =  0; 

zg    =  linspace (0.001, 0.025, 500); 

xgb  =  xg-xb; 
zgb  =  zg-zb; 
theta  =  atan(-xgb./zgb) ; 

alpha  =  m-zwdot; 
beta  =  -zw; 
Al  =  m-xudot; 
Bl  =  -2*cd; 
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A  =  (m-zwdot)*(iy-mqdot)-(m*xg+zqdot)*(mwdot+m*xg) ; 
B   =  -zw*(iy-mqdot)-(m-zwdot)*(mq-m*xg)-(zq+m)*(mvdot+m*xg) . . . 
- (m*xg+zqdot ) *mw ; 

for  i  =  l:length(zg) 

deltaz(i)  =  xgb*w*sin(theta(i) )-zgb(i)*w*cos(theta(i)) ; 

C(i)   =  -deltaz(i)*(m-zvdot)+zv*(mq-m*xg)-(zq+m)*mw; 

D(i)   =  deltaz(i)*zw; 

A2(i)  =  (m*zg(i))~2*(m-zwdot); 

B2(i)  =  -(m*zg(i))~2*zv; 

CA  =  -A1*A; 

CB  =  -A1*B+B1*A; 

CC  =  -Al*C(i)+Bl*B; 

CD  =  -Al*D(i)+Bl*C(i); 

CE  =  Bl*D(i); 

p  =  [CA  CB  CC  CD  CE]  ; 
r(l:4,i)  =  roots(p) ; 

lamda(l:4,i)  =  r(l :4,i)-(r(l :4,i) . ~3.* (alpha. *r(l :4,i)+beta)) 
.*(m.*zg(i))."2./(4.*CA.*r(l:4,i)."3+. . . 
3.*CB.*r(l:4,i).~2+2.*CC.*r(l:4,i)+CD); 

PA  =  -A1*A+A2(i); 
PB  =  -A1*B+B1*A+B2(i); 
PC  =  -Al*C(i)+Bl*B; 
PD  =  -Al*D(i)+Bl*C(i); 
PE  =  Bl*D(i); 

p4  =  [PA  PB  PC  PD  PE]  ;     '/,  e- values  of  4  by  4  system. 
r4(l:4,i)  =  roots(p4) ; 
end 
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7,  Program  crit.O.m 

7.  Evaluation  of  critical  speed  for  delta=0 


rho  = 

1. 

94; 

g   ■ 

32 

..2; 

L 

12 

.9792; 

ndl  = 

0. 

5*rho*L~2 

nd2  = 

0. 

5*rho*L~3 

nd3  = 

0. 

5*rho*L"4 

nd4  = 

0. 

5*rho*L~5 

iy  = 

561.32/nd4; 

zqdot 

= 

-6.33e-4; 

zwdot 

s 

-1.4529e-2; 

zq 

= 

7.545e-3; 

zv 

= 

-1.391e-2; 

mqdot 

= 

-8.8e-4; 

mwdot 

= 

-5.61e-4; 

mq 

= 

-3.702e-3; 

mw 

= 

1.0324e-2; 

cd 

= 

0.015; 

zb 

= 

0/L; 

m 

= 

1556.2363/(g*nd2); 

xudot 

= 

-0.05*m; 

xb 

= 

0/L; 

xg  =  linspace(-0.01,0.01,40) ; 
zg  =  linspace(0. 005, 0.025, 40) ; 

Gv  =  1  -  nm.*(zq+m) ./(zw.*(mq-m.*xg)) ; 
lambda  =  -2*cd/m-xudot ;    7.  lambda  =  -1.6439 

xgb  =  xg-xb; 
zgb  =  zg-zb; 

for  j = 1 : length ( zg ) 
for  i=l:length(xg) 

theta(i,j)  =  atan(-xgb(i)/zgb( j)) ; 

aO  =   (m-zwdot)*(iy-mqdot)-(mwdot+m*xg(i))*(zqdot+m*xg(i)) 
bO  =   (-zwdot*m-m*mv-zq*m)*xg(i)+(-m*mq+zwdot*mq-zqdot*mw. 
-zq*mwdot-m*mwdot-iy*zw+mqdot*zw) ; 
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cO  =  -m*zw*xg(i)+mq*zw-zq*mw-m*mw; 

cl  =   (-m*xg(i)+zwdot*xg(i)+m*xb-zwdot*xb)*6in(theta(i,  j)) . 
+(-m*zb-zwdot*zg( j ) +zwdot*zb+m*zg( j ) ) *cos (theta( i , j ) ) ; 
dl  =   (zw*xg(i)-zw*xb)*sin(theta(i, j))+. . . 
(zw*zb-zw*zg(j))*cos(theta(i, j)) ; 

'/.  After  applying  AD=BC... 

w(i,j)  =  bO*cO/(aO*dl-bO*cl); 

uO(i,j)     =  (1556.2363/(ndl*v(i,j)))-.5; 

end 
end 
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7,  File  crit_delta.m 

7.  Evaluation  of  critical  speed  for  nonzero  stern  plane  angle. 


rho  = 

1. 

94; 

g   = 

3: 

!.2; 

L 

13 

1.9792; 

ndl  = 

0. 

5*rho*L~2 

nd2  = 

0. 

5*rho*L~3 

nd3  = 

0. 

5*rho*L~4 

nd4  = 

0. 

5*rho*L~5 

iy  = 

561.32/nd4; 

zqdot 

= 

-6.33e-4; 

zwdot 

= 

-1.4529e-2; 

zq 

= 

7.545e-3; 

zv 

= 

-1.391e-2; 

mqdot 

= 

-8.8e-4; 

mvdot 

= 

-5.61e-4; 

mq 

= 

-3.702e-3; 

mv 

= 

1.0324e-2; 

m 

= 

1556.2363/(g*nd2); 

xudot 

= 

-0.05*m; 

xb 

= 

0/L; 

zb 

= 

0/L; 

zg 

= 

0.015; 

cdz 

= 

0.5;             7. 

if  c 

zd 

= 

-5.603e-03; 

md 

= 

-2.409e-03; 

EO 

= 

0.1036509;       7. 

E0  = 

El 

= 

-2.3629982e-03;    7. 

El  = 

E2 

= 

7.367214" 

re-03 ;    7. 

E2  = 

7.  if  cdz=0,  division  by  zero  provides  NaN 


=  Integral  of  b(x)*dx 
=  Integral  of  x*b(x)*dx 
=  Integral  of  x~2*b(x)*dx 


xg  =  linspace(-0.01,0.01,80) ; 
uO  =  linspace(2,5,40) ; 


7.  NON-DIMENSIONAL 


xgb  =  xg-xb; 
zgb  =  zg-zb; 
wl    =  1556. 2363./ (ndl. *u0, 


•2); 


delta  =  input ('Enter  the  value  of  delta:  delta  =  linspace( 
delta  =  delta*pi/180; 


)'); 
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for  k=l :length(delta) 
for  j=l:longth(uO) 
for  i=l: length (xg) 

if  delta(k)  >  0 
wO  =  -(zw+sqrt(zw~2-4*cdz*E0*zd*delta(k)))/(2*cdz*E0) ; 
coeffl  =  (mw*wO-cdz*wO'*2*El+md*delta(k))/wl(j); 

polyl   =  [(xgb(i)~2+zgb~2)  (-2*coeff l*zgb)  (coeff l~2-xgb(i)~2)] ; 
thetaO  =  asin(roots(polyl)) ; 

checkl  ■  xgb(i)*cos(thetaO(l))+zgb*sin(thetaO(l)) ; 
check2  =  xgb(i)*cos(theta0(2) )+zgb*sin(th©ta0(2)) ; 
if  abs(check2-coeff 1)  <  abs(checkl-cooff 1) 

thetaO  =  theta0(2) ; 
else 

thetaO  =  thetaO(l); 
end 

B=[  (m-zwdot)     -(m*xg(i)+zqdot)  0 
-(mvdot+m*xg(i))   iy-mqdot     0 
0  0  1]; 

A=[(zv+2*cdz*E0*w0)      (zq+m)-2*cdz*El*wO         0 
mv-2*cdz*El*wO   (mq-m*xg(i))-m*zg*wO+2*cdz*E2*wO  ... 

(xgb(i)*sin(thetaO)-zgb*cos( thetaO ))*wl(j) 
0  1  0]; 

elseif  delta(k)  <  0 
wO  =  -(zv+sqrt(zw~2+4*cdz*E0*zd*delta(k)))/(-2*cdz*E0); 
coeffl  =  (mv*wO+cdz*wO~2*El+md*delta(k))/wl(j); 

polyl   =  [(xgb(i)"2+zgb"2)  (-2*coeff l*zgb)  (coeff l~2-xgb(i)~2)] ; 
thetaO  =  asin(roots(polyl)) ; 

checkl  =  xgb(i)*cos(thetaO(l))+zgb*sin(thetaO(l)); 
check2  =  xgb(i)*cos(theta0(2))+zgb*sin(theta0(2)) ; 
if  abs(check2-coeff 1)  <  abs(checkl-coeff 1) 

thetaO  =  thetaO (2) ; 
else 

thetaO  =  thetaO(l) ; 
end 

B=[  (m-zwdot)     -(m*xg(i)+zqdot)  0 
-(mwdot+m*xg(i))   iy-mqdot     0 
0  0  1]; 
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A=[(zw-2*cdz*E0*w0)       (zq+m)+2*cdz*El*wO         0 
mw+2*cdz*El*wO   (raq-m*xg(i))-m*zg*wO-2*cdz*E2*wO  ... 
(xgb(i)*sin(thetaO)-zgb*cos(thetaO))*wl(j) 
0  1  0]; 

end 

evals  =  eig(A,B) ; 

degstab(i.j)  =  max (real (evals)) ; 

end 
end 

7,  Evaluation  of  the  critical  speed,  ucrit: 

for  i=l : length (xg) 
for  j=l:(length(uO))-l 

flagl  =  sign(degstab(if j)) ; 

flag2  =  sign(degstab(i, j+1)) ; 

if  flagl  '=  flag2 

ucrit(i,k)  =  (uO(j+l)*degstab(i,j)-uO(j)*degstab(i, j+1))/ . . . 
(degstab(i, j)-degstab(i, j+1)) ; 

end 
end 
end 
end 

7,  Revaluation  of  the  nominal  angle,  'thetacr' 

7.  at  the  critical  speed,  'ucrit': 

7.  Note  that  in  this  case,  at  a  certain  'xg' , 

7.  corresponds  a  specific  'ucrit'. 

7,  That's  why  we  use  the  same  index  'i': 

w2  =  1556. 2363. /(ndl .*ucrit . ~2) ; 
for  i= 1 : length ( xg ) 

if  delta(l)  >  0 
wO  =  -(zw+sqrt(zw~2-4*cdz*E0*zd*delta(l)))/(2*cdz*E0) ; 
coeff2  =  (mw*wO-cdz*wO~2*El+md*delta(l))/w2(i) ; 

poly2   =  [(xgb(i)~2+zgb"2)  (-2*coeff 2*zgb)  (coeff2~2-xgb(i)~2)] ; 
thetaO  =  asin(roots(poly2)) ; 
check!  =  xgb(i)*cos(thetaO(l))+zgb*sin(thetaO(l)) ; 
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check2  =  xgb(i)*cos(theta0(2))+zgb*sin(theta0(2) ) ; 
if  abs(check2-coeff2)  <  abs(checkl-coeff 2) 

thetacr(i)  =  thetaO(2); 
else 

thetacr(i)  =  thetaO(l); 
end 
elseif  delta(l)  <  0 
wO  =  -(zv+sqrt(zw~2+4*cdz*E0*zd*delta(l)))/(-2*cdz*E0) ; 
coeff2  =  (mw*wO+cdz*wO~2*El+md*delta(l))/w2(i) ; 

poly2   =  [(xgb(i)~2+zgb~2)  (-2*coeff 2*zgb)  (coeff 2"2-xgb(i)"2)] ; 
thetaO  =  asin(roots(poly2)) ; 

checkl  =  xgb(i)*cos(thetaO(l))+zgb*sin(thetaO(l)) ; 
check2  =  xgb(i)*cos(theta0(2) )+zgb*sin(theta0(2) ) ; 
if  abs(check2-coeff2)  <  abs(checkl-coeff2) 

thetacr(i)  =  theta0(2); 
else 

thetacr(i)  =  thetaO(l); 
end 
end 

end 

results  =  [xg'  ucrit(:,l)  wO*ones(length(xg) ,1)  thetacr'] 
save  ucrit0.dat  results  -ascii 
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C     PROGRAM  HOPF.O.FOR 

C     EVALUATION  OF  HOPF  BIFURCATION  FORMULAS 

C     ZERO  DIVE  PLANE  ANGLE 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  L,IY,MASS,MQDOT,MWDOT,NDl .THETA, 

1  MQ,MW,MD,MDS,MDB,K1,K2, 

2  ALFA, BETA, GAMA, DELTA, 

3  E0,E1,E2,E3,E4, 

4  DW1,DW2,DW3,DW4, 

5  DQ1,DQ2,DQ3,DQ4 

DOUBLE  PRECISION  Mil ,M12,M13,M21 .M22.M23 , 

1  M31.M32.M33, 

2  N11,N12,N13,N21,N22,N23, 

3  N31.N32.N33, 

4  L21 , L22 , L23 , L24 , L31 , L32 , L33 , L34 , 

5  L25,L26,L27,L35,L36,L37, 

6  L21A,L22A,L23A,L24A,L31A, 

7  L32A.L33A.L34A 
C 

DIMENSION  A(3,3) ,T(3,3) ,TINV(3.3) ,FV1(3) ,IV1(3) ,YYY(3,3) 
DIMENSION  WR(3),WI(3),TSAVE(3,3),TLUD(3,3),IVLUD(3),SVLUD(3) 
DIMENSION  ASAVE(3,3) ,A1(3,3) ,A2(3,3) ,XL(25) ,BR(25) 
DIMENSION  VEC0(25) ,VEC1(2S) ,VEC2(25) ,VEC3(25) ,VEC4(25) 

C 

OPEN  (20 , FILE= ' HOPF . RES ' , STATUS= ' NEW ' ) 

C 

WEIGHT=15S6.2363 


IY 

=   561.32 

L 

=      13.9792 

RHO 

1.94 

WRITE 

(*,*)    '    ENTER  CD' 

READ 

(*,*)    CD 

CD 

=  0.5*CD*RH0 

G 

=      32.2 

XB 

0.0 

WRITE 

(*,*)    '    ENTER  ZG/L 

READ 

(*,*)    ZG 

ZG 

=ZG*L 

ZB 

0.0 

MASS 

=  WEIGHT/G 
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BOY 

WEIGHT 

ND1 

0.5*RH0*L**2 

ZqDOT  =■ 

-6.3300E-04+0 

5*RH0*L**4 

ZWDOT  =■ 

■1.4529E-02*0 

5*RH0*L**3 

zq 

7.5450E-03+0 

5*RH0*L**3 

zw 

-1.3910E-02+0 

S*RH0*L**2 

MqDOT  =■ 

-8.8000E-04*0 

5*RH0*L**5 

MWDOT  =■ 

-5.6100E-04*0 

5*RH0*L**4 

Mq 

-3.7020E-03*0 

.5*RH0*L**4 

MW 

1.0324E-02*0 

.5*RH0*L**3 

c 

ZGB=ZG-ZB 

c 

c 

DEFINE  THE  LENGTH  X  l 

IND  BREADTH 

c 

XL(  1)  = 

0.0000 

XL(  2)  = 

0.1000 

XL(  3)  = 

0.2000 

XL(  4)  = 

0.3000 

XL(  5)  = 

0.4000 

XL(  6)  = 

0.5000 

XL(  7)  = 

0.6000 

XL(  8)  = 

0.7000 

XL(  9)  = 

1.0000 

XL(10)= 

2.0000 

XL(11)= 

3.0000 

XL(12)= 

4.0000 

XL(13)= 

7.7143 

XL(14)= 

10.0000 

XL(15)= 

15.1429 

XL(16)= 

16.0000 

XL(17)= 

17.0000 

XL(18)= 

18.0000 

XL(19)= 

19.0000 

XL(20)= 

20.0000 

XL(21)= 

20.1000 

XL(22)= 

20.2000 

XL(23)= 

20.3000 

XL(24)= 

20.4000 

XL(25)= 

20.4167 

DO    102   N   =    1,25 
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102 


XL(N)    =    (L/20.)*XL(N)    -   L/2, 
CONTINUE 


BR(  1)=  0 

000 

BR(  2)=  0 

485 

BR(  3)=  0 

658 

BR(  4)=  0 

778 

BR(  5)=  0 

871 

BR(  6)=  0 

945 

BR(  7)=  1 

010 

BR(  8)=  1 

060 

BR(  9)=  1 

180 

BR(10)=  1 

410 

BR(11)=  1 

570 

BR(12)=  1 

660 

BR(13)=  1 

670 

BR(14)=  1 

670 

BR(15)=  1 

670 

BR(16)=  1 

630 

BR(17)=  1 

370 

BR(18)=  0 

919 

BR(19)=  0 

448 

BR(20)=  0 

195 

BR(21)=  0 

188 

BR(22)=  0 

168 

BR(23)=  0 

132 

BR(24)=  0 

053 

BR(25)=  0 

000 

DO  104  K  ■ 

=  1,25 

VECO(K)= 

=BR(K) 

VEC1(K)= 

=XL(K)*BR(K) 

VEC2(K)= 

=XL(K)*XL(K)*BR(K) 

VEC3(K)= 

=XL(K)*XL(K)*XL(K)*BR(K) 

VEC4(K)= 

=XL(K)*XL(K)*XL(K)*XL(K)*BR(K) 

104  CONTINUE 

CALL  TRAP(25.VEC0 

,XL,EO) 

CALL  TRAP(25,VEC1 

,XL,E1) 

CALL  TRAP(25,VEC2 

,XL,E2) 

CALL  TRAP(25,VEC3 

,XL,E3) 

CALL  TRAP( 

:25,VEC4 

,XL,E4) 
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WRITE  (*,*)  '  ENTER  GAMMA1 

READ   (*,*)  EPSILON 

XGMIN=-0.01 

XGMAX=+0.01 

IXG=80 

XGMIN=XGMIN*L 

XGMAX=XGMAX*L 


DO  1  IT=1,IXG 

WRITE  (*,3001)  IT.IXG 

XG=XGMIN+(XGMAX-XGMIN)*(IT-1)/(IXG-1) 

XGB=XG-XB 

DV=(MASS-ZWDOT) + (IY-MqDOT) - (MASS+XG+ZQDOT) * (MASS+XG+MWDOT) 

CD6=CD/(3.DO*EPSIL0N*DV) 

DW1=CD6*( (IY-MQDOT) *(-EO)+(MASS*XG+ZQDOT)*El) 

DW2=CD6*( (IY-MqDOT) *(3*E1)-(MASS*XG+ZQD0T)*3*E2) 

DW3=CD6*( (IY-MQDOT) *(-3*E2)+(MASS*XG+ZqD0T)*3*E3) 

DW4=CD6*( (IY-MqDOT) *(E3)-(MASS*XG+ZqD0T)*E4) 

Dqi=CD6*((MASS-ZWD0T)*(El)+(MASS*XG+MWD0T)*(-E0)) 

DQ2=CD6*((MASS-ZWD0T)*(-3*E2)+(MASS*XG+MWD0T)*(3*E1)) 

Dq3=CD6*((MASS-ZWD0T)*(3*E3)+(MASS*XG+MWD0T)*(-3*E2)) 

Dq4=CD6*((MASS-ZWDOT)*(-E4)+(MASS*XG+MWDOT)*(E3)) 

THETAO=ATAN(-XGB/ZGB) 

AAO=(MASS-ZWDOT)*(IY-MqDOT)-(MWDOT+MASS*XG)*(ZqDOT+MASS*XG) 

BBO=(-ZWDOT*MASS-MASS*MW-Zq*MASS)*XG 
ft  +(-MASS*Mq+ZWDOT*Mq-ZqDOT*MW 

ft  -ZQ*MWDOT-MASS*MWDOT-IY*ZW+MqDOT*ZW) 

CCO=-MASS*ZW*XG+Mq*ZW-Zq*MW-MASS*MW 

CC1=(-MASS*XG+ZWDOT*XG+MASS*XB-ZWDOT*XB)*SIN(THETAO) 
ft  +(-MASS*ZB-ZWDOT*ZG+ZWDOT*ZB+MASS*ZG)*COS(THETAO) 

DDl=(ZW*XG-ZW*XB)*SIN(THETAO)+(ZW*ZB-ZW*ZG) 
ft  *COS(THETAO) 

C 

WEI=BBO*CCO/(AAO*DDl-BBO*CCl) 

U0=DSqRT(1556.2363/WEI) 

U=UO 
C 

C       DETERMINE  [A]  AND  [B]  COEFFICIENTS 
C 

AllDV=(IY-MqDOT)*ZW+(MASS*XG+ZQDOT)*MW 

A12DV= ( IY-MqDOT) * (MASS+Zq) + (MASS*XG+ZqDOT) * (Mq-MASS*XG) 

A13DV=- (MASS*XG+ZqDOT) *WEIGHT 
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A21DV=(MASS-ZWD0T)*MW+(MASS*XG+MWD0T)*ZW 

A22DV=(MASS-ZWD0T)*(MQ-MASS*XG)+(MASS*XG+MWD0T)*(MASS+Zq) 

A23DV=- (MASS-ZWDOT) +WEIGHT 
C 

A11=A11DV/DV 

A12=A12DV/DV 

A13=A13DV/DV 

A21=A21DV/DV 

A22=A22DV/DV 

A23=A23DV/DV 
C 

C11DV=(IY-MQD0T)*MASS*ZG 

C12DV=- (MASS+XG+ZQDOT) *MASS*ZG 

C21DV=- (MASS-ZWDOT) *MASS*ZG 

C22DV=(MASS*XG+MWDOT) *MASS*ZG 
C 

C11=C11DV/DV 

C12=C12DV/DV 

C21=C21DV/DV 

C22=C22DV/DV 

C       EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 
C 

Kl=-(XGB*SIN(THETAO)-ZGB*COS(THETAO)) 

K2=-(1./6.)*(ZGB*COS(THETAO)-XGB*SIN(THETAO)) 
C 

A(1,1)=0.0 

A(l,2)=0.0 

A(l,3)=1.0 

A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22*U 

DO  11  1=1,3 
DO  12  J=l,3 

ASAVE(I,J)=A(I,J) 
12     CONTINUE 
11    CONTINUE 

CALL  RG(3,3,A,WR,WI,1,YYY,IV1,FV1,IERR) 

CALL  DSOMEG(IEV,WR,WI, OMEGA, CHECK) 


68 


0MEGAO=OMEGA 
DO  5  1=1,3 

T(I,1)=  YYY(I.IEV) 

T(I,2)=-YYY(I,IEV+1) 

5  CONTINUE 

IF  (IEV.EQ.l)  GO  TO  13 
IF  (IEV.Eq.2)  GO  TO  14 
STOP  3004 
14    DO  6  1=1,3 

T(I,3)=YYY(I,1) 

6  CONTINUE 
GO  TO  17 

13    DO  16  1=1,3 

T(I,3)=YYY(I,3) 

16  CONTINUE 

17  CONTINUE 
C 

C       NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 
C 

IN0RM=1 

IF  (INORM.NE.O)  CALL  NORMAL (T) 
C 

C       INVERT  TRANSFORMATION  MATRIX 
C 

DO  2  1=1,3 
DO  3  J=l,3 

TINV(I,J)=0.0 
TSAVE(I,J)=T(I,J) 

3  CONTINUE 
2    CONTINUE 

CALL  DLUD(3,3,TSAVE,3,TLUD,IVLUD) 
DO  4  1=1,3 

IF  (IVLUD(I).EQ.O)  STOP  3003 

4  CONTINUE 

CALL  DILU(3,3,TLUD,IVLUD,SVLUD) 
DO  8  1=1,3 
DO  9  J=l,3 

TINV(I,J)=TLUD(I,J) 
9     CONTINUE 
8    CONTINUE 
C 
C       CHECK  Inv(T)*A*T 
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IMULT=1 

IF  (IMULT.Eq.l)  CALL  MULT(TINV,ASAVE,T,A2) 

IF  (IMULT.EQ.O)  STOP 

P=A2(3,3) 

PEIG=P 
C 

C       DEFINITION  OF  Nij 
C 

N11=TINV(1,1) 

N21=TINV(2,1) 

N31=TINV(3,1) 

N12=TINV(1,2) 

N22=TINV(2,2) 

N32=TINV(3,2) 

N13=TINV(1,3) 

N23=TINV(2,3) 

N33=TINV(3,3) 
C 

C       DEFINITION  OF  Mij 
C 

M11=T(1.1) 

M21=T(2,1) 

M31=T(3,1) 

M12=T(1,2) 

M22=T(2,2) 

M32=T(3,2) 

M13=T(1,3) 

M23=T(2,3) 

M33=T(3,3) 
C 

C       DEFINITION  OF  Lij 
C 

L2E=C11*M31*M31+C12*M21*M31 

L26=2*C11*M31*M32+C12*(M21*M32+M22*M31) 

L27=C11*M32*M32+C12*M22*M33 

L35=C22*M31*M31+C21*M21*M31 

L36=2*C22*M31*M32+C21*(M21*M32+M22*M31) 

L37=C22*M32*M32+C21*M33*M22 
C 

C       DEFINITION  OF  ALFA,  BETA,  GAMA 
C 
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Dl  =N32*L25  +  N33*L3S 
D2  =N32*L26  +  N33+L36 
D3  =N32*L27  +  N33+L37 

Dll=-P 

D12=0MEGA0 

D21=-2*0MEGA0 

D22=-P 

D23=2*0MEGA0 

D32=-0MEGA0 

D33=-P 

BETA=(D2-D21*D1/D11-D23*D3/D33)/(D22-D21*D12/D11-D23*D32/D33) 

ALFA=(D1-D12*BETA)/D11 

GAMA=(D3-D32*BETA)/D33 

L21A=2*C11*ALFA*M31*M33+C12*ALFA*(M21*M33+M23*M31) 

L22A=2*C11*ALFA*M32*M33    +   2*C11*BETA*M31*M33 
+      C12*ALFA*(M22*M33+M32*M23) 
+     C12*BETA*(M21*M33+M23*M31) 

L23A=2*C11*GAMA*M31*M33+2*C11+BETA*M32*M33 
+      C12*GAMA*(M21*M33+M23*M31) 
+      C12*BETA*(M22*M33+M23*M32) 

L24A=2*C11*GAMA*M32*M33+C12*GAMA*(M22*M33+M23*M32) 

L31A=2*C22*ALFA*M31*M33+C21*ALFA*(M21*M33+M23*M31) 

L32A=2*C22*ALFA*M32*M33+2*C22*BETA*M31*M33 
+      C21*ALFA*(M22*M33+M32*M23) 
+      C21*BETA*(M21*M33+M23*M31) 

L33A=2*C22*GAMA*M31*M33+2*C22*BETA*M32*M33 
+      C21*GAMA*(M21*M33+M23*M31) 
+      C21*BETA*(M22*M33+M23*M32) 

L34A=2*C22*GAMA*M32*M33+C21*GAMA*(M22*M33+M23*M32) 

L21=L21A+A13*K2*M11**3+DW1*M21**3+DW2*M31*M21**2 
+  DW3*M21*M31**2+DW4*M31**3 
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L22=L22A+3*A13*K2*M12*M11**2+3*DW1*M22*M21**2 
ft     +     DW2*(2*M21*M22*M31+M32*M21**2) 
ft     +     DW3*(2*M21*M31*M32+M22*M31**2) 
ft     +     3*DW4*M32*M31**2 

L23=L23A+3*A13*K2*M11*M12**2+3*DW1*M21*M22**2 
ft     +     DW2*(M31*M22**2+2*M21*M22*M32) 
ft     +     DW3*(M21*M32**2+2*M22*M31*M32) 
ft     +     3*DW4*M31*M32**2 

L24=L24A+A13*K2*M12**3+DW1*M22**3+DW2*M32*M22**2 
ft     +     DW3*M22*M32**2+DW4*M32**3 

L31=L31A+A23*K2*Mll**3+Dqi*M21**3+Dq2*M31*M21**2 
k  +     Dq3*M21*M31**2+DQ4*M31**3 

L32=L32A+3*A23*K2*M12*Mll**2+3*Dqi*M22*M21**2 
ft     +     Dq2*(2*M21*M22*M31+M32*M21**2) 
ft     +     Dq3*(2*M21*M31*M32+M22*M31**2) 
ft     +     3*Dq4*M32*M31**2 

L33=L33A+3*A23*K2*Mll*M12**2+3*Dqi*M21*M22**2 
ft     +     Dq2*(M31*M22**2+2*M21*M22*M32) 
ft     +     Dq3*(M21*M32**2+2*M22*M31*M32) 
ft     +     3*Dq4*M31*M32**2 

L34=L34A+A23*K2*M12**3+Dqi*M22**3+Dq2*M32*M22**2 
ft     +     Dq3*M22*M32**2+Dq4*M32**3 
C 

R11=N12*L21+N13*L31 

R12=N12*L22+N13*L32 

R13=N12*L23+N13*L33 

R14=N12*L24+N13*L34 

R21=N22*L21+N23*L31 

R22=N22*L22+N23*L32 

R23=N22*L23+N23*L33 

R24=N22*L24+N23*L34 
C 

C       EVALUATE  DALPHA  AND  DOMEGA 
C 

UINC=0.001 

UR  =U+UINC 

UL  =U-UINC 

U  =UR 
C 

A(1,1)=0.0 

A(l,2)=0.0 

A(l,3)=1.0 
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A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22*U 

CALL  RG(3,3,A,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
ALPHR=DEOS 
OMEGR=FREq 

U=UL 

A(1,1)=0.0 

A(l,2)=0.0 

A(l,3)=1.0 

A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22*U 

CALL  RG(3,3,A,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREQ) 
ALPHL=DEOS 
OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) / (UR-UL) 

DOMEGA= ( OMEGR- OMEGL ) / (UR-UL ) 
C 

C       EVALUATION  OF  HOPF  BIFURCATION  COEFFICIENTS 
C 

C0EF1=3 . 0+R11+R13+R22+3 . 0+R24 

C0EF2=3.0*R21+R23-R12-3.0*R14 

AMU2  =-C0EFl/(8.0*DALPHA) 

BETA2=0.25*COEF1 
C       TAU2  =-(C0EF2-D0MEGA*C0EFl/DALPHA)/(8.O*0MEGAO) 
C       PER   =2.0*3. 1415927/OMEGAO 
C       PER   =PER*U/L 

WRITE  (20,2001)  XG/L,U,COEFi ,DALPHA,0MEGAO ,PEIG 
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1   CONTINUE 
STOP 

1001  FORMAT  ('    ENTER   NUMBER   OF   DATA   LINES') 

1002  FORMAT  ('    ENTER  UO ,    ZG,    AND   DSAT') 

1003  FORMAT  ('    ENTER  BOW  PLANE  TO    STERN   PLANE   RATIO') 

1004  FORMAT  ('    ENTER   ZG') 
2001    FORMAT  (6E14.5) 
4001    FORMAT  (3F15.5) 
3001   FORMAT  (215) 

END 

"<—  —  —  —  —  —  —  —  —  —  —  —  —  =  —  _  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —  —.  —  —  —  —  —  —  —  —  —  ... 

SUBROUTINE   DSOMEG  (  UK ,  WR ,  WI ,  OMEGA ,  CHECK  ) 
IMPLICIT   DOUBLE   PRECISION    (A-H.O-Z) 
DIMENSION   WR(3),WI(3) 
CHECK=-1.0D+25 
DO    1   1=1,3 

IF    (WR(I).LT. CHECK)    GO   TO    1 

CHECK=WR(I) 

IJ=I 
1   CONTINUE 

OMEGA=DABS(WI(IJ)) 

IF    (WI(U).GT.O.DO)    IJK=IJ 

IF    (WI(IJ).LT.O.DO)    IJK=IJ-1 

RETURN 

END 

^___=  =  =  =  — _  =  _  — —  —  —  —  —  — =  =  — —  =  =  — —  —  =  =  =  =  — =  =  =  =  =  =  — _=  =  ==  =  =  =  — =  -=  =  =  =  =  . 

SUBROUTINE   DSTABL(DEOS,WR,WI .OMEGA) 
IMPLICIT   DOUBLE   PRECISION    (A-H.O-Z) 
DIMENSION   WR(3),WI(3) 
DE0S=-1.0D+20 
DO    1   1=1,3 

IF    (WR(I).LT.DEOS)    GO   TO    1 

DEOS=WR(I) 

IJ=I 
1   CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS( OMEGA) 
RETURN 
END 

*'—    —    —   —  —  —  —  ————  —  —  —  —  —  —  H2.--  =  =  _  =  --<.  =  =*  =  ---  =  =  ==--  =  -  =  --<.  =  -  =  =  l 

SUBROUTINE   NORMAL (T) 

IMPLICIT  DOUBLE  PRECISION    (A-H.O-Z) 
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DIMENSION  T(3,3),TN0R(3,3) 
CFAC=T(1,1)**2+T(1,2)**2 
IF  (DABS(CFAC).LE.(l.D-lO))  STOP  4001 
TN0R(1,1)=1.D0 

TN0R(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 
TN0R(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 
TN0R(l,2)=O.DO 

TN0R(2.2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 
TN0R(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
DO  1  1=1,3 
DO  2  J=l,2 

T(I,J)=TNOR(I,J) 
2   CONTINUE 

1  CONTINUE 
RETURN 
END 

SUBROUTINE  MULT(TINV,A,T,A2) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  TINV(3,3) ,A(3,3) ,T(3,3) ,A1(3,3) ,A2(3,3) 
DO  1  1=1,3 
DO  2  J=l,3 

A1(I,J)=0.D0 

A2(I,J)=0.D0 

2  CONTINUE 
1  CONTINUE 

DO  3  1=1,3 
DO  4  J=l,3 
DO  5  K=l,3 

A1(I,J)=A(I,K)*T(K,J)+A1(I,J) 

5  CONTINUE 
4   CONTINUE 

3  CONTINUE 
DO  6  1=1,3 

DO  7  J=l,3 
DO  8  K=l,3 

A2(I,J)=TINV(I,K)*A1(K,J)+A2(I,J) 
8     CONTINUE 
7    CONTINUE 

6  CONTINUE 

DO  11  1=1,3 
C        WRITE  (*,101)  (A(I,J),J=1,3) 
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11  CONTINUE 

DO  12  1=1,3 
C        WRITE  (*,101)  (T(I,J),J=1,3) 

12  CONTINUE 

DO  10  1=1,3 
C       WRITE  0,101)  (A2(I,J),J=1,3) 

10  CONTINUE 
C      WRITE  (*,101)  A2(l,l) 
RETURN 
101  FORMAT  (3E15.5) 
END 

SUBROUTINE  TRAP(N,A,B,OUT) 
C 

C     NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

DIMENSION  A(1),B(1) 

N1=N-1 

OUT=0 . 0 

DO  1  1=1, Nl 

OUT1=0.5*(A(I)+A(I+1))*(B(I+1)-B(I)) 
OUT  =0UT+0UT1 
1  CONTINUE 
RETURN 
END 
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C     PROGRAM  HOPF.DELTA.FOR 

C     EVALUATION  OF  HOPF  BIFURCATION  FORMULAS 

C     STERN  PLANE  ANGLE  IS  NON  ZERO 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DOUBLE  PRECISION  L,IY,MASS ,MqDOT,MWDOT,NDl ,THETA, 

1  MQ,MW,MD,ZD,MDS,MDB,K1,K2, 

2  ALFA, BETA, GAMA, DELTA, 

3  E0,E1,E2,E3,E4, 

4  DW1,DW2,DW3,DW4, 

5  DQl,DQ2,Dq3,DQ4 

DOUBLE  PRECISION  Mil .M12.M13 ,M21 .M22.M23 , 

1  M31.M32.M33, 

2  N11,N12,N13,N21,N22,N23, 

3  N31,N32,N33, 

4  L21,L22,L23,L24,L31,L32,L33,L34, 

5  L25,L26,L27,L35,L36,L37, 

6  L21A,L22A,L23A,L24A,L31A, 

7  L32A.L33A.L34A 
C 

DIMENSION  A(3,3) ,T(3,3) ,TINV(3,3) ,FV1(3) ,IV1(3) ,YYY(3,3) 
DIMENSION  WR(3) ,WI(3) ,TSAVE(3,3) ,TLUD(3,3) ,IVLUD(3) ,SVLUD(3) 
DIMENSION  ASAVE(3,3),A1(3,3),A2(3,3),XL(25),BR(25) 
DIMENSION  VEC0(25) ,VEC1(25) ,VEC2(2S) ,VEC3(2S) ,VEC4(25) 

C 

OPEN  ( 10 , FILE= ' HOPF . DAT ' , STATUS= ' OLD ' ) 
OPEN  (20, FILE=' HOPF. RES' , STATUS= ' NEW ) 

C 

WEIGHT=1556.2363 


IY 

=   561.32 

L 

=      13.9792 

RHO 

1.94 

WRITE 

(*,*)    '    ENTER  CD' 

READ 

(*,*)    CD 

CD 

=   0.5*CD*RH0 

G 

=      32.2 

XB 

0.0 

WRITE 

(*,*)    '    ENTER  ZG/L 

READ 

(*,*)    ZG 

ZG 

=ZG*L 

ZB 

0.0 
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MASS   = 

WEIGHT/G 

BOY 

WEIGHT 

ND1 

0.5*RH0*L**2 

ZqDOT  =- 

■6.3300E-04*0 

5*RH0*L**4 

ZWDOT  =■ 

-1.4529E-02+0 

5*RH0*L**3 

ZQ 

7.5450E-03+0 

5*RH0*L**3 

ZW 

■1.3910E-02+0 

5*RH0*L**2 

MqDOT  =- 

■8.8000E-04*0 

,5*RH0*L**5 

MWDOT  =- 

-5.6100E-04*0 

.5*RH0*L**4 

Mq 

-3.7020E-03*0 

.5*RH0*L**4 

c 

MW 

1.0324E-02+0 

. 5*RH0*L**3 

ZGB=ZG-ZB 

c 

c 
c 

DEFINE  THE  LENGTH  X  j 

yjD  BREADTH 

XL(  1)  = 

0.0000 

XL(  2)  = 

0.1000 

XL(  3)  = 

0.2000 

XL(  4)  = 

0.3000 

XL(  5)  = 

0.4000 

XL(  6)  = 

0.5000 

XL(  7)  = 

0.6000 

XL(  8)= 

0.7000 

XL(  9)  = 

1.0000 

XL(10)= 

2.0000 

XL(11)= 

3.0000 

XL(12)= 

4.0000 

XL(13)= 

7.7143 

XL(14)= 

10.0000 

XL(1S)= 

15.1429 

XL(16)= 

16.0000 

XL(17)= 

17.0000 

XL(18)= 

18.0000 

XL(19)= 

19.0000 

XL(20)= 

20.0000 

XL(21)= 

20.1000 

XL(22)= 

20 . 2000 

XL(23)= 

20.3000 

XL(24)= 

20.4000 

XL(25)= 

20.4167 
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DO  102  N  =  1,25 

XL(N)  =  (L/20.)*XL(N)  -  L/2, 
CONTINUE 


BR( 

1) 

BR( 

2) 

BR( 

3) 

BR( 

4) 

BR( 

5) 

BR( 

6) 

BR( 

7) 

BR( 

8) 

BR( 

9) 

BR( 

10) 

BR( 

11) 

BR( 

12) 

BR( 

13) 

BR( 

!14) 

BR( 

:i5) 

BR( 

:i6) 

BR( 

:i7) 

BR( 

:i8) 

BR( 

19) 

BR( 

:20) 

BR( 

.21) 

BR( 

22) 

BR( 

23) 

BR( 

24) 

BR( 

.25) 

.000 
.485 
.658 
.778 
.871 
.945 
.010 
.060 
.180 
.410 
.570 
.660 
.670 
.670 
.670 
.630 
.370 
.919 
.448 
.195 
.188 
.168 
.132 
.053 
.000 


DO  104  K 

VECO(K) 

VECl(K) 

VEC2(K) 

VEC3(K) 

VEC4(K) 

104  CONTINUE 

CALL  TRAP 

CALL  TRAP 

CALL  TRAP 

CALL  TRAP 


=  1,25 

=BR(K) 

=XL(K)*BR(K) 

=XL(K)*XL(K)*BR(K) 

=XL(K)*XL(K)*XL(K)*BR(K) 

=XL(K)*XL(K)*XL(K)*XL(K)*BR(K) 

(25,VEC0,XL,E0) 
(25,VEC1,XL,E1) 
(25,VEC2,XL,E2) 
(25,VEC3,XL,E3) 
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c 


CALL  TRAP(25,VEC4,XL,E4) 
WRITE  (*,*)  '  ENTER  GAMMA1 
READ   (*,*)  EPSILON 
IXG=80 


DO  1  IT=1,IXG 

READ  (10,*)  XG,U,WO,THETAO 

WRITE  (*,3001)  IT.IXG 

XG=XG*L 

XGB=XG-XB 

DV=(MASS-ZWDOT)*(IY-MQDOT)-(MASS*XG+ZQDOT)*(MASS*XG+MWDOT) 

C0NST2=CD/ (3 . 0*DV*EPSILON) 

TDWl=C0NST2*((IY-MqD0T)*(-EO)+(MASS*XG+ZqD0T)*(El)) 

TDW2=CONST2*((IY-MqDOT)*(3.0*El)+(MASS*XG+ZQDOT)*(-3.0*E2)) 

TDW3=C0NST2*((IY-MqD0T)*(-3.O*E2)+(MASS*XG+ZqD0T)*(3.O*E3)) 

TDW4=C0NST2*((IY-MqD0T)*(E3)+(MASS*XG+ZqD0T)*(-E4)) 

TDqi=C0NST2*((MASS-ZWD0T)*(El)+(MASS*XG+MWD0T)*(-EO)) 

TDq2=CONST2*((MASS-ZWDOT)*(-3.0*E2)+(MASS*XG+MWDOT)*(3.0*El)) 

TDq3=CONST2*((MASS-ZWDOT)*(3.0*E3)+(MASS*XG+MWDOT)*(-3.0*E2)) 

TDq4=C0NST2* ( (MASS-ZWDOT) * ( -E4) + (MASS+XG+MWDOT) * (E3 ) ) 

C 

C0NST=TANH(WO/EPSIL0N)*CD/DV 

DW1=C0NST*(  (IY-MqD0T)*(-EO)     +(MASS*XG+ZqDOT)*El) 
DW2=C0NST*(  (IY-MqD0T)*(2*El)    -(MASS*XG+ZqD0T)+2*E2) 
DW3=C0NST*(  (IY-MqD0T)*(-E2)     +(MASS*XG+ZqD0T)*E3) 
Dqi=CONST*(  (MASS-ZWD0T)*(E1)    +(MASS*XG+MWD0T)*(-EO)) 
Dq2=C0NST*(  (MASS-ZWD0T)*(-2*E2)+(MASS*XG+MWD0T)*2*E1) 
Dq3=C0NST*(  (MASS-ZWDOT) *(E3)    +(MASS*XG+MWD0T)*(-E2) ) 

C 

C  DETERMINE    [A]    AND    [B]    COEFFICIENTS 


AllDV=(IY-MqD0T)*(ZW-2*CD*EO*WO) 

+(MASS*XG+ZqD0T)*(MW+2*CD*El*WO) 
A12DV=(IY-MqD0T)*(MASS+Zq+2*CD*El*WO)+(MASS*XG+ZqD0T) 

*(Mq-MASS*XG-MASS*ZG*W0-2*CD*E2*W0) 
A13DV=- (MASS*XG+ZqDOT) *WEIGHT 
A21DV=(MASS-ZWD0T)*(MW+2*CD*E1*W0) 

+(MASS*XG+MWD0T)*(ZW-2*CD*E0*WO) 
A22DV=(MASS-ZWD0T)*(Mq-MASS*XG-MASS*ZG*WO-2*CD*E2*WO) 

+(MASS*XG+MWDOT)*(MASS+Zq+2*CD*E1*WO) 
A23DV=- (MASS-ZWDOT) +WEIGHT 
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A11=A11DV/DV 

A12=A12DV/DV 

A13=A13DV/DV 

A21=A21DV/DV 

A22=A22DV/DV 

A23=A23DV/DV 
C 

C11DV=(IY-MQD0T)*MASS*ZG 

C12DV=-(MASS*XG+ZqD0T)*MASS*ZG 

C21DV=- (MASS-ZWDOT) *MASS*ZG 

C22DV=(MASS*XG+MWD0T)*MASS*ZG 
C 

C11=C11DV/DV 

C12=C12DV/DV 

C21=C21DV/DV 

C22=C22DV/DV 

C       EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 
C 

K1=-(XGB*SIN(THETAO)-ZGB*COS(THETAO)) 
K2=-(1./6.)*(ZGB*COS(THETAO)-XGB*SIN(THETAO)) 
C 

A(1,1)=0.0 
A(l,2)=0.0 
A(l,3)=1.0 
A(2,1)=A13*K1 
A(2,2)=A11*U 
A(2,3)=A12*U 
A(3,1)=A23*K1 
A(3,2)=A21*U 
A(3,3)=A22*U 
DO  11  1=1,3 
DO  12  J=l,3 

ASAVE(I,J)=A(I,J) 
12     CONTINUE 
11    CONTINUE 

CALL  RG(3,3,A,WR,WI,1,YYY,IV1,FV1,IERR) 
CALL  DSOMEG( IEV , WR , WI , OMEGA , CHECK) 
OMEGAO=OMEGA 
DO  5  1=1,3 

T(I,1)=  YYY(I.IEV) 
T(I,2)=-YYY(I,IEV+1) 


81 


5  CONTINUE 

IF  (IEY.EQ.l)  GO  TO  13 
IF  (IEV.EQ.2)  GO  TO  14 
STOP  3004 
14    DO  6  1=1,3 

T(I,3)=YYY(I,1) 

6  CONTINUE 
GO  TO  17 

13    DO  16  1=1,3 

T(I,3)=YYY(I,3) 

16  CONTINUE 

17  CONTINUE 
C 

C       NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 
C 

IN0RM=1 

IF  (INORM.NE.O)  CALL  NORMAL (T) 
C 

C       INVERT  TRANSFORMATION  MATRIX 
C 

DO  2  1=1,3 
DO  3  J=l,3 

TINV(I,J)=0.0 
TSAVE(I,J)=T(I,J) 

3  CONTINUE 
2    CONTINUE 

CALL  DLUD(3,3,TSAVE,3,TLUD,IVLUD) 
DO  4  1=1,3 

IF  (IVLUD(I).EQ.O)  STOP  3003 

4  CONTINUE 

CALL  DILU(3,3,TLUD,IVLUD,SVLUD) 
DO  8  1=1,3 
DO  9  J-1,3 

TINV(I,J)=TLUD(I,J) 
9     CONTINUE 
8    CONTINUE 
C 

C  CHECK   Inv(T)*A*T 

C 

IMULT=1 

IF    (IMULT.EQ.l)    CALL   MULT (TINV, ASA VE.T.A2) 

IF    (IMULT.EQ.O)    STOP 
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P=A2(3,3) 

PEIG=P 
C 

C  DEFINITION   OF   Nij 

C 

N11=TINV(1,1) 

N21=TINV(2,1) 

N31=TINV(3,1) 

N12=TINV(1,2) 

N22=TINV(2,2) 

N32=TINV(3,2) 

N13=TINV(1,3) 

N23=TINV(2,3) 

N33=TINV(3,3) 
C 

C       DEFINITION  OF  Mij 
C 

M11=T(1,1) 

M21=T(2,1) 

M31=T(3,1) 

M12=T(1,2) 

M22=T(2,2) 

M32=T(3,2) 

M13=T(1,3) 

M23=T(2,3) 

M33=T(3,3) 
C 

C       DEFINITION  OF  Lij 
C 

L25=C11*M31*M31+C12*M21*M31 

L26=2*C11*M31*M32+C12*(M21*M32+M22*M31) 

L27=C11*M32*M32+C12*M22*M33 

L35=C22*M31*M31+C21*M21*M31 

L36=2*C22*M31*M32+C21* (M21*M32+M22*M31 ) 

L37=C22*M32*M32+C21*M33*M22 
C 

C       DEFINITION  OF  ALFA,  BETA,  GAMA 
C 

Dl  =N32*L25  +  N33+L35 

D2  =N32*L26  +  N33+L36 

D3  =N32*L27  +  N33+L37 
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Dll=-P 

D12=0MEGA0 

D21=-2*0MEGAO 

D22=-P 

D23=2*0MEGA0 

D32=-0MEGA0 

D33=-P 

BETA=(D2-D21*Dl'r"    D23*D3/D33)/(D22-D21*D12/D11-D23*D32/D33) 

ALFA=(D1-D12*BETA)/D11 

GAMA=(D3-D32*BETA)/D33 

L21A=2*(C11+DW3)^LFA*M31*M33+(C12+DW2)*ALFA*(M21*M33+M23*M31) 
ft  +2*DW1*M21*M23*ALFA 

L22A=2*(C11+DW3)*ALFA*M32*M33    +   2*(C11+DW3)*BETA*M31*M33 
ft  +      (C12+DW2)*ALFA*(M22*M33+M32*M23) 

k  +      (C12+DW2)*BETA*(M21*M33+M23*M31) 

k  +      2*DW1*(M21*M23*BETA+M22*M23*ALFA) 

L23A=2*(C11+DW3)*GAMA*M31*M33+2*(C11+DW3)*BETA*M32*M33 
k  +      (C12+DW2)*GAMA*(M21*M33+M23*M31) 

ft  +      (C12+DW2)*BETA*(M22*M33+M23*M32) 

ft  +      2*DW1*(M21*M23*GAMA+M22*M23*BETA) 

L24A=2*(C11+DW3)*GAMA*M32*M33+(C12+DW2)*GAMA*(M22*M33+M23*M32) 
ft  +2*DW1*M22*M23+GAMA 

L31A=2*(C22+DQ3)*ALFA*M31*M33+(C21+DQ2)*ALFA*(M21*M33+M23*M31) 
ft  +2*Dqi*M21*M23*ALFA 

L32A=2*(C22+DQ3)*ALFA*M32*M33    +   2*(C22+DQ3)*BETA*M31*M33 
ft  +      (C21+DQ2)*ALFA*(M22*M33+M32*M23) 

ft  +      (C21+DQ2)*BETA*(M21*M33+M23*M31) 

ft  +      2*Dqi*(M21*M23*BETA+M22*M23*ALFA) 

L33A=2*(C22+DQ3)*GAMA*M31*M33+2*(C22+DQ3)*BETA*M32*M33 
ft  +      (C21+DQ2)*GAMA*(M21*M33+M23*M31) 

ft  +      (C21+DQ2)*BETA*(M22*M33+M23*M32) 

ft  +      2*DQ1*(M21*M23*GAMA+M22*M23*BETA) 

L34A=2*(C22+DQ3)*GAMA*M32*M33+(C21+DQ2)*GAMA*(M22*M33+M23*M32) 
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ft  +2*Dqi*M22*M23*GAMA 

C 

L21=L21A+A13*K2*M11**3 

L22=L22A+3*A13*K2*M12*M11**2 

L23=L23A+3*A13*K2*M11*M12**2 

L24=L24A+A13*K2*M12**3 

L31=L31A+A23*K2*M11**3 

L32=L32A+3*A23*K2*M12*M11**2 

L33=L33A+3*A23*K2*M11*M12**2 

L34=L34A+A23*K2*M12**3 
C 

L21=L21+TDW1*M21**3+TDW2*M31*M21**2+TDW3*M21*M31**2+TDW4*M31**3 

L22=L22+TDW1*3.0*M22*M21**2+TDW2+(2.0*M21*M31*M22+M32*M21**2) 
ft  +TDW3*(2.0*M31*M32*M21+M22*M31**2)+TDW4*3.0*M32*M31**2 

L23=L23+TDW1*3.0*M21*M22**2+TDW2*(M31*M22**2+2.0*M21*M22*M32) 
ft  +TDW3* (M21*M32**2+2 . 0*M31*M32*M22) +TDW4+3 . 0*M31*M32**2 

L24=L24+TDW1*M22**3+TDW2*M32*M22**2+TDW3*M22*M32**2+TDW4*M32**3 

L31=L31+TDQl+M21**3+TDq2*M31*M21**2+TDQ3*M21*M31**2+TDQ4*M31**3 

L32=L32+TDqi*3.0*M22*M21**2+TDq2*(2.0*M21*M31*M22+M32*M21**2) 
ft  +TDq3*(2.0*M31*M32*M21+M22*M31**2)+TDq4*3.0*M32*M31**2 

L33=L33+TDqi*3.0*M21*M22**2+TDq2*(M31*M22**2+2.0*M21*M22*M32) 
ft  +TDq3*(M21*M32**2+2.0*M31*M32*M22)+TDq4*3.0*M31*M32**2 

L34=L34+TDqi*M22**3+TDq2*M32*M22**2+TDq3*M22*M32**2+TDq4*M32**3 
C 

R11=N12*L21+N13*L31 

R12=N12*L22+N13*L32 

R13=N12*L23+N13*L33 

R14=N12*L24+N13*L34 

R21=N22*L21+N23*L31 

R22=N22*L22+N23*L32 

R23=N22*L23+N23*L33 

R24=N22*L24+N23*L34 
C 

C       EVALUATE  DALPHA  AND  DOMEGA 
C 

UINC=0.001 

UR  =U+UINC 

UL  =U-UINC 

U  =UR 
C 

A(1,1)=0.0 

A(l,2)=0.0 
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A(l,3)=1.0 

A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22*U 

CALL  RG(3,3,A,WR,WI,0,YYY,IV1,FV1,IERR) 
CALL  DSTABL(DEOS,WR,WI,FREq) 
ALPHR=DEOS 
OMEGR=FREQ 

U=UL 

A(1,1)=0.0 

A(l,2)=0.0 

A(l,3)=1.0 

A(2,1)=A13*K1 

A(2,2)=A11*U 

A(2,3)=A12*U 

A(3,1)=A23*K1 

A(3,2)=A21*U 

A(3,3)=A22*U 
C 

CALL  RG(3,3,A,WR,WI,0,YYY,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHL=DEOS 

OMEGL=FREq 
C 

DALPHA=(ALPHR-ALPHL)/(UR-UL) 

DOMEGA= ( OMEGR- OMEGL ) / (UR-UL ) 
C 

C       EVALUATION  OF  HOPF  BIFURCATION  COEFFICIENTS 
C 

C0EF1=3 . 0*Rll+R13+R22+3 . 0+R24 

C0EF2=3 . 0*R21+R23-R12-3 . 0*R14 

AMU2  =-C0EFl/(8.0*DALPHA) 

BETA2=0.25*COEF1 
C       TAU2  =-(C0EF2-D0MEGA*C0EFl/DALPHA)/(8.0*0MEGA0) 
C       PER   =2.0*3. 1415927/OMEGAO 
C       PER   =PER*U/L 
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WRITE  (20,2001)  XG/L ,U,C0EF1 ,DALPHA,0MEGAO ,PEIG 
1  CONTINUE 
STOP 

1001  FORMAT  ('  ENTER  NUMBER  OF  DATA  LINES') 

1002  FORMAT  ('  ENTER  UO,  ZG,  AND  DSAT') 

1003  FORMAT  ('  ENTER  BOW  PLANE  TO  STERN  PLANE  RATIO') 

1004  FORMAT  ('  ENTER  ZG') 
2001  FORMAT  (6E14.5) 
4001  FORMAT  (3F15.5) 
3001  FORMAT  (215) 

END 

SUBROUTINE  DSOMEG ( UK , WR ,  WI ,  OMEGA ,  CHECK ) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  WR(3),WI(3) 
CHECK=-1.0D+25 
DO  1  1=1,3 

IF  (WR(I).LT. CHECK)  GO  TO  1 

CHECK=WR(I) 

IJ=I 
1  CONTINUE 

OMEGA=DABS(WI(IJ)) 
IF  (WI(IJ).GT.O.DO)  IJK=IJ 
IF  (WI(IJ).LT.O.DO)  IJK=IJ-1 
RETURN 
END 
0================================================================= 

SUBROUTINE  DSTABL ( DEOS , WR , WI , OMEGA ) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  WR(3),WI(3) 
DEOS=-1.0D+20 
DO  1  1=1,3 

IF  (WR(I).LT.DEOS)  GO  TO  1 

DEOS=WR(I) 

IJ=I 
1  CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS( OMEGA) 
RETURN 
END 

SUBROUTINE  NORMAL (T) 
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IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  T(3,3),TN0R(3,3) 
CFAC=T(1,1)**2+T(1,2)**2 
IF  (DABS(CFAC).LE.(l.D-lO))  STOP  4001 
TN0R(1,1)=1.D0 

TN0R(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 
TN0R(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 
TN0R(l,2)=O.DO 

TN0R(2,2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 
TN0R(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
DO  1  1=1,3 
DO  2  J=l,2 

T(I,J)=TNOR(I,J) 
2   CONTINUE 

1  CONTINUE 
RETURN 
END 

SUBROUTINE  MULT(TINV,A,T,A2) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 
DIMENSION  TINV(3,3),A(3,3),T(3,3),A1(3,3),A2(3,3) 
DO  1  1=1,3 
DO  2  J=l,3 

A1(I,J)=0.D0 

A2(I,J)=0.D0 

2  CONTINUE 
1  CONTINUE 

DO  3  1=1,3 
DO  4  J=l,3 
DO  5  K=l,3 

A1(I,J)=A(I,K)*T(K,J)+A1(I,J) 

5  CONTINUE 
4   CONTINUE 

3  CONTINUE 
DO  6  1=1,3 

DO  7  J=l,3 
DO  8  K-1,3 

A2(I,J)=TINV(I,K)*A1(K,J)+A2(I,J) 
8     CONTINUE 
7   CONTINUE 

6  CONTINUE 

DO  11  1=1,3 
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C        WRITE  (*,101)  (A(I,J),J=1,3) 

11  CONTINUE 

DO  12  1=1,3 
C        WRITE  (+.101)  (T(I,J),J=1,3) 

12  CONTINUE 

DO  10  1=1,3 
C       WRITE  (*,101)  (A2(I,J),J=1,3) 

10  CONTINUE 
C      WRITE  (*,101)  A2(l,l) 
RETURN 
101  FORMAT  (3E15.5) 
END 


C= 


SUBROUTINE  TRAP(N ,A ,B,OUT) 
C 

C     NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

DIMENSION  A(1),B(1) 

N1=N-1 

0UT=O . 0 

DO  1  1=1, Nl 

OUT1=0.5*(A(I)+A(I+1))*(B(I+1)-B(I)) 
OUT  =0UT+0UT1 
1  CONTINUE 
RETURN 
END 
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C     PROGRAM  SIM. FOR 

C 

C     SIMULATION  PROGRAM  USING  A  FOURTH  ORDER  ADAMS-BASHFORTH 

C     INTEGRATION  SCHEME 

C 

REAL  L,MW,MASS,IY,MqDOT,MWDOT,MQ, 

*  U,TIME,ND1,ND2,ND3,ND4, 

*  F1(4),F2(4),F3(4) 


DIMENSION  XL(25) ,BR(2S) ,VEC1(25) ,VEC2(2S) 


OPEN  (10,FILE='SIM.DAT' ,STATUS='OLD' ) 

OPEN  (20,FILE='SIM.RES' ,STATUS='NEW  ) 
C 

C     ENTER  SPEEDS  AND  TIME  DATA 
C 

READ(10,*)  U,TSIM,DELT,IPRNT, 
*  ZG,XG,THETA,W,Q 

C 

C     GEOMETRIC  PROPERTIES  AND  HYDRODYNAMIC  COEFFICIENTS 
C 

PI     =4.0*ATAN(1.0) 

THETA  =THETA*PI/180 

RHO        =1.94 

CDZ        =0 . 50 

L  =13.9792 

ND1        =0.5*RHO*L**2 

ND2        =O.5*RH0*L**3 

ND3        =0.5*RHO*L**4 

ND4        =0.5*RHO*L**5 

WEIGHT=1556 . 2363/ (ND1*U**2) 

MASS   =1556. 2363/ (32. 2+ND2) 

IY    =561.32/ND4 

ZQDOT=-6.3300E-04 

ZWD0T=-1.4529E-02 

ZQ    =  7.5450E-03 

ZW    =-1.3910E-02 

MQDOT=-8.8000E-04 

MWDOT=-5.6100E-04 

Mq    =-3.7020E-03 
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MW    =  1.0324E-02 

c 

ZB     =0.0 

XB     =0.0 

ZGB    =ZG-ZB 

XGB    =XG-XB 

c 

c 

DETERMINE  [A]  AND  [B]  COEFFICIENTS 

c 

DV    =(MASS-ZWDOT) * (IY-MqDOT) - (MASS+XG+ZQDOT) * (MASS+XG+MWDOT) 

A11DV=(IY-MQD0T)*ZW+(MASS*XG+ZQD0T)*MW 

A12DV=(IY-MqD0T)*(MASS+ZQ)+(MASS*XG+ZQD0T)*(MQ-MASS*XG) 

A13DV=- (MASS+XG+ZQDOT) +WEIGHT 

A21DV=(MASS-ZWD0T)*MW+(MASS*XG+MWD0T)*ZW 

A22DV= (MASS-ZWDOT) * (MQ-MASS+XG) + (MASS*XG+MWDOT ) * (MASS+Zq ) 

A23DV=- (MASS-ZWDOT) *WEIGHT 


A11=A11DV/DV 

A12=A12DV/DV 

A13=A13DV/DV 

A21=A21DV/DV 

A22=A22DV/DV 

A23=A23DV/DV 

c 

c 

DEFINE  THE  LENG' 

c 

XL(  1)  = 

0.0000 

XL(  2)  = 

0.1000 

XL(  3)  = 

0.2000 

XL(  4)  = 

0.3000 

XL(  5)  = 

0.4000 

XL(  6)  = 

0.5000 

XL(  7)  = 

0.6000 

XL(  8)  = 

0.7000 

XL(  9)= 

1.0000 

XL(10)= 

2.0000 

XL(11)  = 

3.0000 

XL(12)= 

4.0000 

XL(13)= 

7.7143 

XL(14)= 

10.0000 

XL(1S)= 

15.1429 

XL(16)= 

16.0000 

X  AND  BREADTH  B  TERMS  FOR  THE  INTEGRATION 
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XL(17)= 
XL(18)= 
XL(19)= 
XL(20)= 
XL(21)= 
XL(22)= 
XL(23)= 
XL(24)= 
XL(25)= 


17.0000 
18.0000 
19.0000 
20.0000 
20.1000 
20.2000 
20.3000 
20.4000 
20.4167 


102 


DO  102  N  =  1,25 

XL(N)  =  (L/20)*XL(N) 
CONTINUE 


-  L/2 


BR(    1) 

= 

0.000 

BR(    2) 

= 

0.485 

BR(    3) 

= 

0.658 

BR(    4) 

■ 

0.778 

BR(    5) 

= 

0.871 

BR(    6) 

= 

0.945 

BR(    7) 

= 

1.010 

BR(   8) 

= 

1.060 

BR(    9) 

= 

1.180 

BR(10) 

= 

1.410 

BR(ll) 

= 

1.570 

BR(12) 

= 

1.660 

BR(13) 

= 

1.670 

BR(14) 

= 

1.670 

BR(15) 

= 

1.670 

BR(16) 

= 

1.630 

BR(17) 

= 

1.370 

BR(18) 

= 

0.919 

BR(19) 

= 

0.448 

BR(20) 

= 

0.195 

BR(21) 

= 

0.188 

BR(22) 

= 

0.168 

BR(23) 

= 

0.132 

BR(24) 

= 

0.053 

BR(25) 

= 

0.000 

DO  103  M  =  1,25 
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XL(M)  =  XL(M)/L 
BR(M)  =  BR(M)/L 
103   CONTINUE 


PISIM  =TSIM/DELT 
ISIM   =PISIM 
C 

C     SIMULATION  BEGINS 
C 

DO  100  1=1, ISIM 
C 

C       CALCULATE  THE  DRAG  FORCE,  INTEGRATE  THE  DRAG  OVER  THE  VEHICLE 
C 

DO  101  K=1,2S 
UCF=W-XL(K)*q 

VEC1(K)=BR(K)*UCF*ABS(UCF) 
VEC2(K)=BR(K)*UCF*ABS(UCF)*XL(K) 
101    CONTINUE 

CALL  TRAP(25,VEC1,XL,HEAVE) 
CALL  TRAP ( 25, VEC2, XL, PITCH) 
HEAVE=CDZ*HEAVE 
PITCH=CDZ*PITCH 
C 

C       COUPLING  EQUATIONS 
C 

DWDV=(IY-MQDOT)*(-HEAVE)+(MASS*XG+ZqDOT)*PITCH 
DQDV=(MASS-ZWDOT)*PITCH+(MASS*XG+MWDOT)*(-HEAVE) 
C1DV=(IY-MQD0T)*MASS*ZG*Q**2-(MASS*XG+ZQD0T)*MASS*ZG*W*Q 
C2DV=-(MASS-ZWD0T)*MASS*ZG*W*Q+(MASS*XG+MWD0T)*MASS*ZG*Q**2 
C 

DW=DWDV/DV 
DQ=DQDV/DV 
C1=C1DV/DV 
C2=C2DV/DV 
C 

WDOT  =A11*W+A12*Q+A13*(XGB*C0S(THETA) 

♦  +ZGB*SIN(THETA))+DW+C1 

qDOT   =A21*W+A22*Q+A23*(XGB*C0S(THETA) 

*  +ZGB*SIN(THETA))+Dq+C2 
THEDOT=q 

IF  (I.GT.3)  GO  TO  150 
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c 

C       INITIAL  FIRST  ORDER  INTEGRATION 
C 

W     =  W     +  DELT+WDOT 

q    =  q    +  DELT*qDOT 

THETA  =  THETA  +  DELT*THEDOT 
C 

Fl(I)=WDOT 

F2(I)=qD0T 

F3(I)=THEDOT 
C 

C       ADAMS -BASHFORTH  INTEGRATION 
C 

150    Fl(4)=All*W+A12*q+A13*(XGB*C0S(THETA)+ZGB*SIN(THETA))+DW+Cl 

F2(4)=A21*W+A22*q+A23*(XGB*COS(THETA)+ZGB*SIN(THETA))+D0+C2 

F3(4)=q 
C 

W    =W    +(DELT/24.0)*(55.0*F1(4)-59.0*F1(3)+37.0*F1(2) 

*  -9.0*F1(1)) 

q    =q     +(DELT/24.0)*(55.0*F2(4)-59.0*F2(3)+37.0*F2(2) 

*  -9.0*F2(1)) 
THETA=THETA+(DELT/24.0)*(55.0*F3(4)-59.0*F3(3)+37.0*F3(2) 

*  -9.0*F3(1)) 

F1(1)=F1(2) 
F1(2)=F1(3) 
F1(3)=F1(4) 
F2(i)=F2(2) 
F2(2)=F2(3) 
F2(3)=F2(4) 
F3(1)=F3(2) 
F3(2)=F3(3) 
F3(3)=F3(4) 

TIME=I*DELT 

J=J+1 

IF  (J.NE.IPRNT)  GO  TO  100 

ALFA=W 

ALFA=(-ATAN(ALFA))*180/PI 

WRITE  (20,991)  TIME, THETA* 180/PI, q , ALFA 

J=0 
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C 


100  CONTINUE 

STOP 
991  FORMAT  (4E15.5) 
END 
C 

SUBROUTINE  TRAP(N,A,B  ,OUT) 
C 

C     NUMERICAL  INTEGRATION  ROUTINE  USING  THE  TRAPEZOIDAL  RULE 
C 

DIMENSION  A(1),B(1) 

N1=N-1 

OUT=0 . 0 

DO  1  1=1, Nl 

0UT1=0.5+(A(I)+A(I+1))*(B(I+1)-B(I)) 
OUT  =0UT+0UT1 
1  CONTINUE 
RETURN 
END 
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